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Abstract 

The study of the properties and dynamics of self-gravitating bosonic objects in 
Einstein gravity was conducted. Bosons are promising candidates for dark matter. 
They can form compact objects through a Jeans instabihty mechanism. We studied 
boson stars made up of self-gravitating scalar fields, with and without nonlinear 
self-couplings. These are non-topological solutions of the coupled Einstein-Klein- 
Gordon equations. We studied the stability of boson stars in the ground and excited 
states, and determined the quasinormal mode (QNM) frequencies of stable boson 
stars in spherical symmetry. The study was carried out in the standard Einstein 
theory of General Relativity and in Brans-Dicke theory. We also studied the for- 
mation of these objects in Brans-Dicke theory showing that they can form from the 
self-gravitation of bosonic matter. We also studied the possibility of a bosonic halo 
surrounding galaxies. These halo models predict the observed flatness of galactic 
curves. We studied their formation and stability. 

After an extensive study in spherical symmetry we carried out numerical studies 
of boson star dynamics in full 3+1 dimension. One focus of the 3 spatial dimension 
(3D) study was on the validation of the numerical code constructed to solve Einstein 
equations with matter sources. The use of the scalar field has unique advantages: 
Boson Stars do not suffer from difficulties associated with hydrodynamic sources 
(like shock waves or the surface problems of neutron stars). They also do not suffer 
from the difficulties related to the singularities of black holes. The code was first 
tested with spherical perturbations and compared with the spherical results. We 
determined the coordinate conditions needed to provide stable evolutions under ra- 
dial perturbations. We then went on to study their behavior under non-spherical 
perturbations. Both scalar and gravitational radiation produced under these per- 
turbations were studied. We reproduced the QNM frequencies of the stars, as 
determined by perturbation studies carried out by other groups. The energy gen- 
erated by the perturbation was studied with different radiation indicators. We also 
observed the collapse to black holes of unstable boson-star configurations. We sim- 
ulated the collision of two boson stars. This is of interest as the two body problem 
is as yet unresolved in general relativity. 
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The total mass of the star is plotted as a function of time. The mass 
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Chapter 1 
Introduction 



The study of bosonic compact objects is motivated by several considerations. From a 
purely aesthetic point of view, the existence of compact objects made up of fermionic 
subparts (like neutron stars), leads one to postulate the existence of bosonic coun- 
terparts made up of particles satisfying Bose-Einstein statistics. This was the mo- 
tivation for the original works of Kaup |jl]], and Ruffini and Bonazzola 0. 

Although not observed thus far, these purely theoretical objects nevertheless can 
provide many useful physical insights. They can be a source of scalar radiation and 
gravitational radiation, and can under certain conditions collapse to black holes. 
Studying the evolution of the resultant black holes then adds a new dimension to 
black hole studies since now matter is present in the spacetime. 

Recently, Bosonic objects have been suggested as strong candidates for dark- 
matter in the universe. This has changed their status from being purely theoretical 
toy models, to possibly being extremely real physical objects. 

They also occupy a significant place in the field of numerical relativity. For the 
first time one can consider the full Einstein system of equations with matter sources 
present. Before the advent of this field, boson stars could only be studied under 
infinitesimal perturbations using perturbative techniques P, Q. 

In one-dimensional numerical relativity we have observed the behavior of boson 
stars under larger and more general spherical perturbations. In 3D, for the first 
time, a fully relativistic problem with matter present was considered. Boson stars 
provide a model which does not suffer from the problematic singularities of black- 
hole spacetimes. In addition the presence of a scalar field with an exponential fall-off 
that, in principle, stretches to spatial infinity, prevents the formation of sharp outer 
surfaces that are characteristic of neutron stars. 

The 3D code could be tested against the well studied ID codes by considering 
spherical perturbations. The gravitational wave signals of non-spherical configura- 
tions in the linear regime were compared against the results of linear non-spherical 
perturbation studies. A comparison of the signals calculated by different methods 
provided a further successful code test. In order to stabilize the code, various nu- 
merical techniques were implemented and studied. With a stable code one could 
then study the signals that came out of large non-spherical perturbations. These 
physical studies were precursors to a planned series of future neutron star studies. 
The collapse of an unstable boson star to a black hole, followed by a "turning on" 



of apparent horizon boundary conditions, was a new aspect of black hole studies 
that boson stars provided. In addition we began simulations of the collision of two 
boson stars. This is of importance in the generalized two-body problem since the 
resulting dynamics may not be sensitively dependent on the constituent objects. 

A study of compact bosonic objects touched upon numerous physical issues. In 
one dimension, the evolution of boson stars in the ground and excited states, with 
and without self-coupling, and the evolution and formation of boson-stars in Brans- 
Dicke theory were investigated. The stability and formation of objects formed from 
massless scalar fields was also studied in a one-dimensional halo model. In 3D the 
study of gravitational waves from the non-spherical perturbation of a boson star was 
followed by a study of black hole evolutions of the black holes formed from collapsed 
boson stars. 

We begin by introducing the Einstein equations with general matter sources, 
followed by a historical perspective on the nature of boson stars and a summary of 
earlier work in the field. The notion of dark matter is then introduced, what makes 
one believe it exists, and why bosons are good candidates for it. This is followed by 
the basics of a Jeans instability mechanism to give an idea of how boson stars could 
actually form. 

The final topics in the introduction are about boson stars without spherical 
symmetry. Since they can be a source of gravitational waves, we first review how 
Einstein's equations predict gravitational radiation in a non-spherically symmetric 
setting. In this regard we also introduce the idea of waveform extraction through 
the Zerilli, Bel-Robinson, and Newman-Penrose functions which we use in our code. 
We discuss some of the other work in the field of Boson stars in higher dimensions 
as well as the results of 3D perturbation studies to which we compare our results. 
The penultimate topic in the introduction is about the 3 + 1 ADM split of space- 
time and how the Einstein equations can be separated into constraint and evolution 
equations. Since our work is numerical, we also present an overview of the kinds of 
difficulties that we face in a general numerical code. The actual details of the code 
and numerics are explained in the chapters that follow. 

1.1 Einstein Equations with Matter: Variational Ap- 
proach 

The central tenet of GR is that matter must respond to geometry by moving in a 
certain way, and geometry must in turn respond to matter by curving. We briefly 
outline the various quantities that we encounter in the equations that we will have 
to solve, and then use a variational approach to derive the equations which a fully 
relativistic system with matter sources must obey H. 



Metric: The invariant measure of proper length or proper time between neigh- 
boring events in spacetime is given by 

ds^ = g^v dx^ dx" = —dr'^, (1) 

where the metric Qfj^^, in a Minkowski spacetime is the diagonal metric r]ij = 
diag{—l, 1, 1, 1). 

F: In curved space variations in tensors are affected by the fact that the basis 
vectors themselves are twisting and turning from spacetime point to spacetime 
point. The Connection Coefficient T measures the changes in the basis vectors 
through the relation 

V,e^ = r-'^, e^. (2) 

Here Vu measures changes along the basis vector 6,^. In any coordinate basis 
([e^,ej,] = ) the connection coefficients are given by 
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t; 9°" i9Sfi,u + 95u,fi - Qixv^i) = r°^^. (3) 
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Geodesic Equation: The path of particles between spacetime points V and Q 
is one that extremizes Jp dr leading to the geodesic equation 

where A is a parameter that parametrizes the path and agrees with the proper 
time r within X = ar + b, with arbitrary constants a and b. 

Riemann tensor: Spacetime curvature results in geodesies deviating from 
nearby geodesies, with test particles accelerating relative to one another. Mea- 
suring the deviation of one geodesic from another involves quantifying this 
curvature through the Riemann curvature tensor. In fiat space (zero curva- 
ture) it must vanish. The Riemann curvature tensor components in terms of 
the connection coefficients are 

dA _ '-^'- a/3 '-^'- afi -pA tct tA -pa /r\ 

"''^ ~ ~d^ dx^ '"'" °^ ~ '"^ "''■ ^^ 

Physically in a curved spacetime parallel transport of a vector around a closed 
loop brings it back to the original point visibly changed. The Riemann tensor 



is a measure of the change. A 4 vector A transported around a closed curve 
formed by vectors u and v (and their commutator to close the gap) suffers a 
change 

(5 A" + R^'f^^s Ai^u"'v^ = 0. (6) 

The Ricci tensor and scalar follow from R^i, = W^^au and R = g^^R^^. 

In contrast to the geometrical quantities above, the stress energy tensor T con- 
tains information about the matter sources of gravity. The conditions of energy and 
momentum conservation are contained in the relation V ■ T = 0. The components 
of T^jy are symmetric. 

The equations that we must rely on to tell us how the system behaves are the 
Einstein equations 

G = 87rT, (7) 

where the geometry of spacetime that tells matter how to move is on the left hand 
while the matter sources telling spacetime how to curve are on the right hand. 
Thus G must only be geometry dependent and so must be entirely made from the 
Riemann tensor and the metric. For simplicity it is made linear in the Riemann 
tensor and since it is an indication of curvature it must vanish in flat space. Like 
the right hand side, it must be a dual tensor that is symmetric. The only second 
rank tensor we can build fulfilling these criteria is (up to an overall constant) 

G^u = R^u — -R Qfiu- (8) 

The factor 8 tt in (7) follows from comparisons to well known situations in Newtonian 
gravity {g^^ ^ Vi^v)- 

It can be shown that we can associate the tt component of the stress-energy tensor 
T*"^ with the density of mass energy, while T^^ is the density of the j component of 
the momentum. The 00 component of the Einstein equation is, therefore, called the 
Hamiltonian constraint and the Oi components are the momentum constraints. 

The Einstein curvature tensor G has to have zero divergence and from this the 
divergence of T^i^ follows. Why this is so can be seen from the fact that a given 
spacetime hypersurface has a metric given by 

ds^ = g^y dx^ dx" . (9) 

Using the metric symmetries we get ten independent components of g^y. When the 
geometry evolves then, if we had ten independent equations 



G^, = 8 7rT^,, (10) 

we could determine all ten g^^, on a slice. This is not possible since we have the 
freedom to choose our coordinates: 

x^-^x'^, (11) 

and this transformation would keep the line element (9) unchanged, but lead to the 
transformation g^^ -^ g' . Therefore four of the ten conditions must be eliminated 
so that we have only six independent Einstein equations. Here is where the vanishing 
divergence of the geometry independent T^,y plays a part. The equation 

T'"'-^ = 0, (12) 

gives us the required four conditions that are subtracted from the ten equations 
leaving six independent geometrical equations. 

In order to actually derive equation (10), a variational principle is applied with 
action written in the form 

1= f Cd^x= f Li-gf^ d^x. (13) 

The Lagrangian density is separated into a geometrical and a field part 

£ = £geom + -Cficld = (-fi')2 L, (14) 

-'-' -'^gcom ~r -infield- 

One then extremizes this action. In the ADM formalism, one finds it convenient to 
think of the dynamics as representing the evolution from the specified 3-geometries 
of an initial spacelike to a final spacelike slice. The action integral is to be extremized 
with regard to choice of spacetime between these faces: 

s^seo, (15) 

with the action changing with 3-geometry according to 

6S = f r^ 5g^j d^x. (16) 

Here tt*-' = -^ is the field momentum conjugate to gij. Consider the simplest action 
principle one can imagine: 

^aeom = [j^j *R- (17) 



Here ^i? is the four-dimensional scalar invariant. To be an invariant, Lgeom niust be 
independent of the choice of coordinate systems. In the neighborhood of a point 
one can always choose coordinates such that the first derivatives of g^^, vanish. 
Therefore, apart from a constant, no scalar invariants can be built homogeneously 
out of the metric coefficients and their first derivatives. Unlike the case of mechanics 
one must, therefore, consider second derivatives of the metric coefficients. One then 
has 20 distinct components of the curvature tensor with six parameters of a local 
Lorentz transformation giving 14 distinct choices. Only, one of these, ^R, is linear 
in second derivatives of metric components. 
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(19) 



Every F should be symmetric in its two lower indices giving a total of 4 x (3*2 + 4) = 
40 independent components. In order for the action to be an extremum, variations 
due to changes in g^^", as well as the Fs' should vanish. 
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5(g^^R^py^) d^x+ 5 
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d X. 



(20) 



The Fs themselves do not transform like tensors although their variations 5T do. 
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dx" dx'^ d'^x'^ 
dx^ dx^ dx^dx^ 



dx^' 



(21) 



The last term destroys the tensor character of the Fs but cancels in 5T^^. From the 
fact that tensors can be studied in any coordinate system, we conveniently choose 
one in which the Fs vanish at the point in question and so from (5) (after replacing 
commas by semicolons for correct formulas in any coordinate system) 
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^r ai3;ij. — ^r Q,^.^, 
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6R, 



afi 



(^r olI3;\ — ^r aA;/3- 



(23) 



Also 



5{-g)^ = o(-^) ^ {-g) 



j^v 
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—^{-9Y9,.h''^- 



(24) 



Thus from (20) we see that 

^ J L"^^ R^fs (-i V^g,. 6g>^''^ + g^^ 5R^p V^ 

+ Rap V-g^a'^^ \<fx + j 5 (Lfieid v^) d'^x, 

= T7 — / ( ^-^ 9aP R + Rap + (^<5r a/3;A — ^^ aA;/3 j 

+ / 5 (L field ^f^) d^x. 
The ST term is evaluated as follows: 

^ I (^"^ (5r^«^;A - 5T\x;p) V^ d^x 

\[g-P5T\p}^-g^P,^5T\p 



(25) 
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(26) 
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Using the fact that the metric is covariantly constant {g°''^-^a = 0) we get 
1 
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"^ (5r\«.x-(5r 



a/3;A 



aA;/3 



'^ d X 



(27) 



I^/({«"^*r^"4.-K*r^-}, 



-g ct'x. 



Using Gauss' theorem we get an integral over the surface. Since variations vanish 
on the boundary this integral identically vanishes leaving 

1 



16 TT 



j y-l^gap R + Rapj V^Sg""^ d^ + J 6 (Lfield V^) d^x = 0. (28) 



The stress-energy tensor is defined to be 

2 I SCf^c\d\ 2 ( 5\/—g -Afield 



TafS = 

Using this definition 

5 ( -Afield 



{-gf^ V ^r^ 



-.al3 



(29) 



/ 5 (Lfieid v^) d'x = I 5g-^^-^^^^^^^ d'x = -\j T^pV^g5g-^ d'x, (30) 



which gives the equations of motion 



Gap = 8nTai3- (31) 

If we now restrict attention to a matter Lagrangian that depends on the metric but 
not metric derivatives (this assumption holds for usual matter fields such as spinors, 
gauge bosons and minimally coupled scalars but fails for non-minimally coupled 
scalars), then 

Gaf3 = 8 Vr I gal3 Lfield — 2 „ ^^ 1 . (32) 

This then sets up the equations for general metrics and matter sources in GR. 
We then specialize them to particular metrics, for example, spherically symmetric 
spacetimes, and a particular matter source that describes bosonic spin zero objects. 
Since most of the work has been on boson stars we start by describing what those 
are. 



1.2 What are Boson Stars? 

Boson stars are compact self-gravitating objects made up of scalars that are held 
together by the balance between the attractive force of gravity and the dispersive 
effects of a wave equation. The scalar field is complex and satisfies the Klein- 
Gordon equation, and the system has field solutions with different numbers of nodes. 
The lowest energy state, or ground state, is characterized by a field distribution 
that has no nodes. Unlike the case of fermions, where one can invoke the perfect 
fluid approximation and write an equation of state, a system of bosons has an 
anisotropic distribution of stress near its center, making such a concept inadmissible. 
Nevertheless, they have many similar properties to neutron stars. For a review, 
seei. 

1.2.1 An Historical Perspective 

From an historical point of view the study of boson stars began with the work of 
Kaup |I|. Electromagnetic geons, as introduced by Wheeler 0] were gravitational- 
electromagnetic objects satisfying the electromagnetic and gravitational field equa- 
tions. Their spherical configurations were found to be in unstable equilibrium. Com- 
pelled by this concept, Kaup introduced a Klein-Gordon geon, which satisfied the 
coupled KG and Einstein equations (the KGE equations) in a spherically symmetric 
metric. From the Lagrangian 

2 
T71 

L = i?+(?'^^0*;^0;,--^0*0, (33) 



with scalar field and a spherically symmetric metric 

dr^ = b^dt^ - r^dn^ - a'dr^ (34) 

is derived a scalar field equation (in appropriately defined units) 

vj/" + (^ + '^' - A') ^' + a^[^-l\-^ = 0, (35) 

where (pir^t) = e*'^*\E'(r), a = C^ and b = e^ . 

Transforming to a new radial coordinate s where 

- dr, 36 

b 

so that ds/dr = a/b, we get 

-^ + — — + (i?^-62)M/ = 0. (37) 

as r a ds ^ ' 

This equation is similar to the radial Schrodinger equation for a potential b^ . In 
the asymptotic region 6^ ~ 1 — 2 Mjr so the potential is similar to the hydrogen 
atom case in that it has a 1/r dependence. This suggests countably infinite localized 
solutions for ^ with different numbers of nodes. The conserved current associated 
with the f/(l) symmetry of the problem is 

J, = ^2(<l>*;,<l> -$;,$*), (38) 

With the particle number times m the mass of a single boson 

Nm = -i l\^ Jo= I r'^abJ^ dr. (39) 

The ground state solutions are then numerically obtained and a plot of mass and 
particle number profiles for different central densities is made. The nature of the 
plot is shown in Fig. 2 (chapter 2). The mass of the bosonic geon increases to a 
maximum and then starts to decrease as does the particle number with its maximum 
at the same value of central density. However Nm is greater than the mass of the 
bosonic geon for central densities up to a point beyond the maximum. 

This result is very interesting since the mass profile as a function of central field 
density is very similar to that of neutron stars and white dwarfs when plotted as a 
function of density p. This is so, as pointed out in this paper, despite the absence of 
an equation of state in the case of a the KGE system where the radial pressure —T^ 
is different from —Tg = —Tf. Not only are the profiles similar, but remarkably, so 
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also are the stability properties, although the usual perfect fluid approach to stellar 
stability can no longer be applied. 

A first-order perturbation approach in spherical symmetry is developed in Kaup's 
paper, whose essential features are that the metric components are perturbed as 
u —^ u + Su and X ^ X + 6X so that 6grr = —2 a^SX where 5z/ and 6X are functions 
of r and t. The scalar field perturbation is of the form 

6<!> = e'^' [R{r,t)+il{r,t)]. (40) 

All time dependences of the perturbations are taken to be of the form e*'^* and 
solutions with imaginary u or negative u'^ are obviously unstable. A variational 
approach follows, to establish the ground state a; and this gives an erroneous result 
that all ground state solutions are stable. 

A Newtonian approach to the problem is presented in 0. The nonrela- 
tivistic Schrodinger equation with a potential V satisfying the Poisson equation 
V'^V = —inGp is used. A mass profile plot as a function of central density 
yields a maximum mass of 1.653 Mpi/m, far greater than the relativistic value of 
0.633 Mpi/m because of the continued use of the Newtonian approach to dense 
configurations no longer in the Newtonian regime. In p[ an A^-body Hamiltonian 
for an N boson configuration interacting via Newtonian gravity is written. The 
hydrogen-like resulting Hamiltonian then yields a ground state expectation value by 
analogy. A lower bound on the ground state energy (better by a factor of almost 
3 than that of [0) is obtained by separating out the center of mass kinetic energy 
(which corresponds to the global motion of the system and should not contribute 
to the ground state energy). This is now much lower than the relativistic result. A 
semi-relativistic calculation with a Hamiltonian of the form 

^=i:(p?+'"^)'-s: ^. (41) 

yields an answer of 1.51 Mpi/m, lower than p| but higher than the relativistic values. 
A Hartree approximation to the problem of A^ identical bosons interacting 
through two-body attractive Yukawa forces with force range - (the gravitational 
case corresponding to ;U = 0) shows bound state solutions for -^ < some critical 
value. In the Hartree approximation the wavefunction of the A^ particle ground 
state is in terms of single particle minimum energy wave functions | / > 

l^0>=|/>l|/>2 |/>7V, (42) 

with the Hamiltonian 
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H = T + V, 



2 m f^^ 



Ev?, 



^ p-fJ-\n-rj\ 



p—h^\'%—']\ 

^ = -9'T. 1 r- (43) 



i>j=i ri '^J 



This paper |^ posits a variational problem in the particle density n{r) = 
N[f*{r)f{r)] with energy minimized as a function of n{r) under the normaliza- 
tion constraint / d^rn{r) = N. For the gravitational case, the results of |^ for 
bound state energy are recovered. 

With the particle number density decaying slowly to infinity, the velocity profile 
of a test particle (calculated from v{r) = JGM{r)/r = JGrn /q n{r') Anr'"^ dr' /r 
by equating centripetal and gravitational forces) as a function of radius does not 
fall off as steeply as the finite boundary fermion case. One of the suggestions of 
this paper was, therefore, to consider boson halos around luminous objects (such as 
galaxies) as a possibility. 

In the literature, there are also treatments of finite temperature boson stars in 
the classical and quantum regimes ||10| , ITTH but our work has involved systems of a 



condensate of zero temperature bosons all of which occupy the same state. 

By far the most thorough analyses on the stability of boson stars from a pertur- 
bative point of view have been by T.D. Lee et. al. While the previous papers have 
dealt only with boson stars in the ground state T.D. Lee et. al. discuss excited state 
configurations as well. In |12[ they discuss the cusp structure of the boson star mass 
and particle number profile in great detail. This cusp structure has important sta- 
bility consequences. The equilibrium configurations have a field with an underlying 
characteristic frequency u so that (p = a{p) e"^*. The action for the system consists 
of a gravitational part and a matter part given by: 

C+oo 

A{g)= dtL{g), (44) 



A{m) = ~j [0*^^0^^ + [/(0*0)] rfr, 

with volume element dr = J\g\dtdpd9 d(j). The metric is taken to be of the form 
ds"^ = —e^'^dt^ + p^dVt^ + e^^dp^ with matter Lagrangian 

L{m) =A7T f {-U -V + W) e"+'' p^ ^^_ (45) 

The total Lagrangian L = /C — V has components 

/C = 4 TT r e"+^ p'W dp = 271 r e^-" p^ ^* ^ dp, (46) 
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where spherical symmetry ensures no kinetic energy is carried off by the gravitational 
field and 



V = V{m) + V{g) = 4 TT / e"+" p^ {U + V) dp- L{g) 

Jo 



(47) 



with [/ = lmV2 + i/V 

2 r 

and y = le-^-^^^' 
2 



dp' 



The Hamiltonian is then 



H = IC + V = IC + V{m) - L{g) = 2 /C - L{m) - L{g) = 2 /C - L. 



(48) 



The particle number N = J j^\g\^dxidx2dx3 is derived from the conserved cur- 



rent discussed earlier in the context of |I[( j° 



dt 



dt 



Lua ) giving 



2K = Nuj and E = Nu~L 



(49) 



and hence 



dM 

In 



dE 

dN 



U!. 



The equilibrium configurations are determined numerically from the Einstein equa- 
tions, with mass determined from the Schwarzschild nature of the asymptotic region: 
namely GM = limp^oo pu- A plot of the mass versus particle number has slope u. 
For any number of nodes n it increases up to a point {M{n, l),N{n, 1)) where it 
has a cusp and abruptly changes course till M(n, 2), N{n, 2) and then another cusp 
causes it to change course again till a third cusp and so on. The sign of dN/du 
changes when a cusp is passed and the function M{N) has a cusp therefore at 



d2 



^ = £^L = 0. At the cusp, in fact, 

dN 
duj 



dM 
dui 



0. 



(50) 



This change of sign at the cusp is an important part of stability discussions in 
where a perturbation technique is used. 

In their stability analysis the perturbed wave function is expressed in terms of 
the radial part of the unperturbed wave function as 



-^e-'''''^[a{p)+XR{p.t)+iXi{p,t)] 



(51) 



where the unperturbed field function satisfies the Klein-Gordon equation written 
for brevity here in operator form ha = 0. The x functions are expanded in terms 
of a complete set of orthonormal radial basis functions fi such that 
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XR{p,t)=J2Ri{t)fi{p) and xi{p,t) = Y.WMp)- (52) 

The perturbed Hamiltonian and Hamilton's equations are written down. The nature 
of the eigenvalues of R and / hold the key to stability. For an n node state it is 
shown that of the eigenvalues A^ of the basis functions fa under operator h, n of them 
are negative using the Sturm-Liouville theorem. For a ground state wave function, 
therefore, all these are positive. The sign of these eigenvalues determines the nature 
of the squared eigenvalues of R and /, being real if A is positive and imaginary if 
not. For a ground state star, therefore, the squared eigenvalues of R and / are real. 
To set about finding out whether the eigenvalues themselves are real (if the squared 
eigenvalues are positive) or imaginary indicating stability or instability respectively, 
a simplifying trick is used. The change of sign only occurs at the cusp hence, if 
any configuration on the first branch is stable, all configurations on that branch 
must be so. Thus a Newtonian configuration on that branch is considered, thereby 
simplifying the equations and it is shown that first branch of ground state boson 
stars are stable. For excited states though, all those configurations are unstable. 
This analysis for stability was done without including any self-coupling terms in the 
matter Lagrangian. 

The need to incorporate the self-coupling term is mostly motivated by astro- 
physical considerations. A single boson star, also called a mini-soliton star, is not 
very large. Although, there are no limits on how many of them there could be it 
would be nice to be dealing with large objects. The radius of the star held up from 
collapse by the Heisenberg uncertainity principle is p ~ 1/R which for a relativistic 
star gives R ~ 1/m where m is the mass of a boson. The mass of the star is of 
order Mpi/m, much smaller than a fermion star of mass of order Mpjm\. To get 
a physical dimension on this, a IGeV boson mass gives us a boson star mass of 
M ~ 10~^^Mq. The importance of the self- interacting term is measured by the ra- 
tio y(0)/m^|0p = A|0p/m^. The energy density of a system without self coupling 
is given by p ~ M/R^ ~ Mpfn?. Since the energy density of non-interacting bosons 
is given by p ~ ?7i^|0p this gives |0| ~ Mpi and hence the self interaction may be 
ignored if A << m? /M'j,^. The mass of the system with self-coupling is of the order 
M ~ yXMch so that even a small self-coupling parameter can significantly increase 
the sizes of these objects as discussed in [ p!3[] . 

In this paper |T^, the equilibrium configurations are set up (with self-coupling 
term) in a spherically symmetric spacetime. The usual mass profile (as a function 
of central field density) for ground-state boson stars with different values of the 
dimensionless self-coupling parameter A = ^2'^^Q is plotted. The maximum mass 
of these profiles steadily increases with increase in A. In the limit of large A it 
is shown that one can write an effective equation of state, and the perfect fluid 
theorems can be used to study their stability, indicating that the branch to the left 
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of the maximum in the mass profile are stable and those to the right are unstable. 
This effective equation of state, and vanishing of anisotropies, makes these high 
A boson stars very much like neutron stars. They have a well defined radius with 
their wave character suppressed. This suggests that one might be able to use the 



self-coupling parameter as an anisotropy control parameter [0]. In realistic systems, 
one expects some anisotropies to be present in systems. These anisotropies could 
affect surface red-shifts and critical masses W^. Usually one puts them in by hand 



in neutron star models in an ad-hoc manner. With self-coupled boson stars this can 
be done more naturally. Interestingly the fractional anisotropy [T^ — Tg)/T^ at the 
"radius" of the star 

is roughly the same for all central densities. 

1.2.2 Boson Stars as Dark Matter Candidates 

In this section we discuss why boson stars, although they have never been seen, 
might indeed really exist. To do so we introduce the notion of dark matter and the 
evidence for this kind of matter. We then show why bosons themselves might be 
good candidates for this kind of matter. We then describe by what mechanism these 
bosons could form compact self-gravitating objects. 

Evidence of Dark Matter 

Almost all the information about the Universe has been obtained through photons 
(radio photons from neutral hydrogen gas. X-ray photons from ionized gas, optical 
photons from stars etc.) ||16|, [l^. What about matter that does not emit detectable 
radiation? For example, the star formation process could result in extremely dim 
main sequence stars with masses below the lower limit for hydrogen burning on the 
main sequence (O.OSM©) which could only be detected at distances of less than Ipc 
with present instruments. 

Even within a given type of astronomical object, there is no reason to expect 
mass and luminosity to be perfectly correlated. Stars fainter than the sun make up 
about 75% of the mass, while stars brighter than the sun account for 95% of the 
luminosity in the solar neighborhood. It is found that most astronomical systems 
have extremely high mass to light ratios. 

Dark matter is defined to be matter whose existence we know of only because 
of gravitational effects. Objects like white dwarfs in the solar neighborhood are not 
dark matter candidates as their existence is inferred from present densities of visible 
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white dwarfs, theories on stellar evolution, and the history of star formation rates 
in the solar neighborhood. 

Zwicky in 1933 was the first to suggest that dark matter existed. His work, based 
on the measurement of radial velocities of 7 galaxies in the Coma cluster, showed 
that individual galaxies had radial velocities differing from the mean velocity of the 
cluster, with an RMS dispersion of 700kms~^ . This dispersion was taken to be a 
measure of the kinetic energy per unit mass for the galaxies in the cluster and a 
crude estimate of the radius of the cluster then provided a measure of the mass 
of the cluster through the virial theorem. The mass to luminosity ratio based on 
this model was found to be almost a factor of 400 times larger than the mass to 
luminosity ratio calculated from the measurement of rotational curves of the nearby 
spirals. He concluded that virtually all the cluster mass was in the form of some 
invisible or dark matter undetectable except by its gravitational force. Although 
based on uncertain cluster radius and distance scales with minimal statistics, it was 
nevertheless a fairly accurate result. The mass to light ratio of the Coma cluster as 
a whole exceeds the mass to light ratio of luminous parts of typical galaxies in the 
cluster by more than a factor of 30. 

Ostriker and Einasto et al. proposed in 1974 that even isolated galaxies had 
large amounts of dark matter around them with spiral galaxies having dark halos 
several times the radius of luminous matter. We have studied one such halo model 
made up of bosonic matter. 

Since boson models have been used to fit rotation curves, and because scalar 
field cosmology is much discussed, we present a few dark matter scenarios in these 
contexts. 

• Galactic Rotation Curves 

From emission lines in HII regions and the 21 cm emission line of neutral hy- 
drogen, the rotation curves of galaxies are optically traced. Using the fact that the 
gravitational force provides the centripetal force 

7 = ^. (54) 

one sees that in the inner region of roughly constant density, where M{r) = p |vrr^, 
the speed v{r) should rise linearly with distance r from the center. An intermediate 
region should exist, where the speed reaches a maximum and starts to decline, until 
the outer Keplerian region where the system acts like a point mass concentrated at 
the center, so v oc r"-*^/^. 

However, most of the rotation curves of the over 70 spiral galaxies studied are 
either fiat or slowly rising up until the last point measured. The few that have falling 
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rotation curves either fall off at a lower than Keplerian rate, or have neighbors that 
may perturb their velocity profiles. This absence of the Keplerian region indicates 
the absence of well determined masses of galaxies, even for those galaxies whose 
rotation curves extend to large enough radii to contain all the light. Thus there 
are no spiral galaxies with accurately determined total mass. This suggests the 
presence of a massive dark halo extending beyond optical radii. A flat rotation 
curve (constant velocity) would indicate halo masses increasing linearly with radius 
out to radii beyond the last observed point. 

• Groups of Galaxies 

Collections of galaxies with separation distances much smaller than typical in- 
tergalactic separations would seem to indicate a gravitational binding between com- 
ponent galaxies. 

If a large number of stars are in a potential 0(x, t) at a given time t the number 
of stars in a given volume d^x with velocities in the range d^v centered on v is given 
by 

/(x,v,t) rf^xrf^v (55) 

where / is the distribution function. The density of stars must satisfy a continuity 
equation since the drift of stars must be star conserving. Hence 

On integrating the above equation over all velocities (after multiplying them by Vj 
and replacing acceleration with the negative gradient of the potential) we get 



Using the divergence theorem, and the fact that no stars are moving infinitely fast, 
the last integral in the above equation when integrated by parts can be modified to 

Also, since the velocity and the coordinate are independent, the partial x deriva- 
tive in the second term can be brought outside the integral. This gives 

^ + ^(^ + ,|i = 0. (59) 

Ot OXi OXj 
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where p = J f(P"v with Vi = - J fvi(P\. Multiplying the above equation by x^ and 
integrating over all space variables, we get 

/ ,,M<i3, + I ,^a(p^ ^3^ ^ I ,,, |i,3, ^ 0. (60) 

In a steady state, the first term vanishes, and one gets the tensor virial theorem 

2K^k + W,k = 0, (61) 

where the kinetic energy tensor K is 

Kjk = -J pvjvl c?^x. (62) 

The kinetic term was obtained by integrating the second integral of (60) by parts. 
The potential energy tensor is given by the third integral in (60). In order to see 
this consider the W tensor 

(63) 
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d^x d^x. (64) 

Differentiating through the second integral (p is a function of the prime coordinate 
and is unaffected by a derivative in the unprimed one), we get 

W,, = G I I p(x) p(x') '"\^'"'~^J^ d^x' d^'x. (65) 

J J |x — x'l"* 

Exchanging the dummy coordinates (x and x') and adding the result to the above 
integral, we get 

W,u = -\g I I P(x) P(x')^^^^^^^^%^^ d^x d^x. (66) 

Taking the trace of both sides of this equation, we get 

W = --G [ p(x) / ^^ rf^x rf^x (67) 

2 J J |x — x'l 

= 2 / P(^) <^(^) ^'^• 
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By considering the change in potential energy of a system, whose density and po- 
tential are p(x) and 0(a;), when we bring in a small mass from spatial infinity to 
position X, we can show that the above expression for W is one of many alterna- 
tive expressions for the potential energy of the system. So we have a scalar virial 
theorem 

2K + W = 0. (68) 

Consider a group of A^ galaxies with masses Mj. Then by the virial theorem 

N N / 1 \ 

^M, <^;2>,= ^^GM,M,(^^) , (69) 

where averages are time averages. Assuming that the mass to light ratios of each 
member galaxy is the same 7 = Mj/Lj and writing the above equation in terms of 
luminosity, we get 



3 71 EliL, <vl^>t,, 



.1=1 J^i ^ ^p,i ^t,n 



3 \|Ri-Rj| 



■j I / t,n 



where we have replaced vf by 3 < v^^ >q, where Vp is the line of sight component 
(which should have the same mean square value as the orthogonal components). 
Also, Ri is the projection of ri onto the plane of the sky, and one can write the 
average rrz^j ^^ 2/7r times the average .^^J;^ , by averaging over the angles. Al- 
though all the temporal and angular quantities cannot be measured, if one has a 
large number A^ of galaxies, then for galaxy orbits of random phases and orientations 
the observable quantity 

_ 3 vr Ell U vl, 

lest ^ ^ j^j 1 ; \'^) 

should approach 7 (and should be a good estimator for 7). If there is dark matter 
present then this estimate would be lower than the actual value. By estimating 
7 using the virial theorem, taking into account a dark matter distribution of the 
same order as the galactic spatial distribution, (note that this would still be smaller 
than the actual matter distribution if there were a more extensive dark matter 
distribution), Huchra and Geller found a median value of 260hMQ in the visible 
band for their groups. The mass to light ratios for groups was seemingly much 
larger than values seen in the luminous parts of galaxies. This suggests again the 
presence of large amounts of dark matter. 
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• Dark Matter Cosmology 

Estimating a lower limit on the density of the universe using the Einstein cosmo- 
logical model, and comparing it to the measured luminosity density, indicates there 
is far more dark matter than visible matter. 

The universe on a large scale appears to be very homogeneous and isotropic, so 
much so that the small scale anisotropies might be considered as perturbations on 
a homogeneous background. In the idealized version, considering total homogeneity 
and isotropy of spatial geometry in Einstein's theory of gravity G^^, = SttT^^, if we 
assume a boundary condition of closure a three sphere would satisfy all conditions. 
The spatial metric of a three sphere can be built up step by step starting from a 1 
sphere. Visualized as embedded in an Euclidean space of one higher dimension, S^ 
satisfies x^ + y^ = a^, which in polar coordinates transforms into x = acos0 and 
y = asin0, giving the metric da"^ = a'^dcjP'. A 2 sphere S*^, x^ + y"^ + z^ = a^, under 
the transformation x = a sin 9 cos (j), y = a sin 6 sin 0, and z = a cos 6 gives a spatial 
metric dcr^ = a'^{d9'^+sin^ Odcp"^). A 3 sphere S^, given by x^+y'^+z'^+w^ = a^, under 
the transformation x = a sin x sin 6 cos (p, y = a sin x sin 6 sin (p, z = a sin x cos 6, and 
w = acosx, therefore, gives da^ = a^i^dx^ + siv? xi^O"^ + sin^^(i0^)). Hence, the 
spacetime geometry is described by 

ds^ = ~df + a^ (dx^ + sin^ x (dO^ + sin^ 9 dcj)'^)) . (72) 

Calculating the connection coefficient symbols F as described above, and the Rie- 
mann tensor from them, one can calculate Ricci tensor and Ricci scalars. The tt 
component of the Einstein equation is then 

3 fda\ 3 „ ,_^, 

Where Tu = p- From red shift measurements and distance measurements one can 
calculate the ratio of the velocity of recession of a galaxy, to the distance to the 
galaxy, which should equal the ratio of the rate of increase in the radius of the 
universe to the radius of the universe. This is the Hubble parameter today: 



1 da 
a dt 



(74) 



From this and (73) 

P>^(|)' = fffo^ (75) 

8 71 a'^ \dt J 8 TT 

From a Hubble time today of -^ (Hubble expansion rate Hq ~ 55km/sec/Mpc) one 
gets a lower limit on the density of pc = 5 x 10~^^ g / cm^ as compared to an observed 
luminosity density of p; ~ 2 x 10~^^ g / cm^ . A large amount of dark matter must be 
present. 
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Bosons as Dark Matter Candidates 

There are suggestions that dark matter is mostly of a non-baryonic nature. About 
10~^ sees after the Big Bang, at a black body temperature of lO^^K, nuclear re- 
actions like n + p ■^^ d + J were in equilibrium. With many high energy photons 
present, the deuterons were as likely to be disassociated as formed, keeping the net 
deuterium density low. Once the temperature cooled to about 10^'' i^, the photon 
energy (kT) was no longer high enough to disassociate the deuterons. The primor- 
dial deuterium densities were "frozen in" and therefore, became a measure of baryon 
density then and now (baryon number conservation). However, whatever deuterons 
formed rapidly burned into He^, H, He^, and Li^; with He^, having the most binding 
energy, dominating. The higher the density of baryons, the faster the deuteron pro- 
ducing reactions; and a higher density of deuterium at the epoch of nucleosynthesis 
would cause more of it to burn off to helium. The observed abundances of deuterium 
(which is really a lower limit on primordial abundances) therefore puts a limit on 
the contribution of baryons to the mass density of the universe. The present mass 
density in baryons calculated comes out to less than 6% of pc- 

As a candidate for structure formation in the universe, bosons are in many 
ways suitable candidates as suggested by |]18|. Spontaneous fluctuations in the pre- 
inflationary epoch could have been greatly magnified by inflation, producing regions 
slightly denser than their surroundings which were amplified by gravity to set up 
the coalescence to present day structures. These fluctuations, or matter density 
variations over the range of clusters of galaxies, can be measured. The square of the 
density fluctuation strength multiplied by the volume over which they are sampled 
provides a power spectrum. If one plots the power spectrum against the length scale 
over which fluctuations are detected for models dominated by cold dark matter and 
hot dark matter, one finds the former is unable to explain large scale structure while 
the latter is unable to explain small scale structure. Hot dark matter candidates 
have low masses and large random velocities. Cold dark matter candidates have large 
masses and low velocities. However, low mass bosons which follow Bose-Einstein 
statistics also have a large number of particles with low velocities unlike neutrinos 
which follow Fermi-Dirac statistics. These are capable of keeping the power peak 
for large scale structures as well as having enough power on small scales. 

Besides this, many particle physics and cosmological models rely on the presence 
of scalars. These scalars have never been seen experimentally and are strong dark 
matter candidates. The actual formation of a boson star out of these scalars must 
rely on a Jeans instability mechanism described in the next subsection. There are 
various dark halo models made up bosonic objects that have been used to fit rotation 
curves. We have investigated the stability and formation of one such model in the 
context of a cosmological model universe that is not closed. 
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Formation of Compact Objects — The Jeans Instability Mechanism 

Although many particle physics and cosmological models predict the existence of 
scalars, the question arises as to how these scalars could come together to form 
compact objects like boson stars. In this regard, we discuss the Jeans instability 
mechanism. From a homogeneous background of matter, local fluctuations can grow 
in time and cause clustering of matter. 

Unlike plasmas that have both positive and negative charges, so as to be neu- 
tral on large scales, an always attractive gravitational force prevents gravitational 
systems from being in static homogeneous equilibrium. In principle an inflnite ho- 
mogeneous gravitational system in equilibrium is impossible. If density and pressure 
are constant and mean velocity vq is zero, then Euler's equation 

9v 1 

+ V.v = — Vp-V0, (76) 

ot p 

leads to V0o = 0. However Poisson's equation gives V^0o = 4:710 po clearly in 
contradiction unless po = 0. In a homogeneous gravitational system, there are 
no pressure gradients to balance the existing gravitational attraction. In order to 
construct an infinite homogeneous gravitational system the "Jeans effect" is used. 

The conditions of the Jeans effect invoke Poisson's equation only when the per- 
turbed density and potential are involved [|^, |^, while the unperturbed potential 
is assumed to be zero (0o = 0; V^0 = 47rG'(p — po))- In uniformly rotating ho- 
mogeneous systems, where one has centrifugal forces in place of pressure gradients 
to balance the equilibrium gravitational field, no Jeans effect is necessary. So a 
homogeneous system can be in static equilibrium in a rotating frame. 

• Physical Basis of the Jeans Instability 

Consider an infinite homogeneous fiuid of density po and pressure po with no internal 
motions so that Vf = 0. Now draw a sphere of radius r around any point and 
compress this region by reducing the volume from V to V — aV. Thus the density 
oc 1/V is perturbed by an amount pi ~ a po, and as a result there is a pressure 
perturbation pi = {dp/dp)opi = v^apQ. The pressure force per unit mass is Fp = 
—Vp/p, and so there is an extra pressure force of magnitude |FpJ = |Vpi/po| ~ 
Pi/(por) ~ av1/r (assuming pi is of the form r" then Vpi goes as pi/r). Also 
because of the increased density there is an extra gravitational force Fci = — V0i. 
Here 0i = G Mjir — 6r) — G M/r ^ G M Sr/r"^. Since 6V = —aV implies 6r ~ —ar 
therefore F^^ ~ G Ma/r"^ with M = ^p^r^. Hence the net force is Fp^ + F^^ 
which if outward, implies the compressed fiuid reexpands and the perturbation is 
stable. If the fiuid continues contracting then the perturbation is unstable with the 
gravitational force larger than the pressure force. This happens if 
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G pqt a > aVg/r, (77) 



or 



v^ 



Perturbations on a scale larger than this are unstable. 

In the case of stellar systems stability can be discussed in the context of this 
mechanism. The density of states / satisfies 

dfi dfi dfo 

^ + v-^-V0,. — = 0. (79) 

Here we have made use of the time independence and homogeneity of the unper- 
turbed density /o(x, v,t) = /o(v) as well as the Jeans effect 0o = 0. Poisson's 
equation gives 

VVi = 4 7rG' I f\d\. (80) 

The solution is /i(x,v, t) = faiy) exp[i(k ■ x— cut)], 0i(x, t) =0aexp[i(k ■ x— cut)] 
where 

(k.v-cc;)/„-0„k- ^ = (81) 

-k^<Pa = ^T^G j fad\. (82) 

Noting that 0i is only a function of x, t and not v and substituting for fa from above 
into the equation for 0^, we get 

r k ■ ^ 
-k^ = 4:nG ^^^d\. (83) 



k ■ \ — uj 

Consider /o to have a Maxwellian velocity distribution 



/o(v) = -^e-^. (84) 



|2 



Taking the x direction to be the direction of k, we get 



^-\-"l/<^'^ 



k a^ j-oo kvx — uj 
where po is the density. For uj = equation (83) becomes 
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dvx = 0. (85) 



e{uj = o) = ej = i^^. (86) 



Instability corresponds to imaginary u, as can be seen from the form of /i. Set 
u = i'-y where 7 is real and positive. After multiplying and dividing by kv^ + i'j, 
(85) becomes 



1,,2 /„2 



' — t^ L -ww^iT """ = »• (87) 

where the term with i 7 in the numerator is odd and hence vanishes (and so has not 
been written). Using J^ x^ exp(— x^) dx/{x^ + /9^) = 0? — n (5 exp(/9^)[l — erf(/9)] 
gives 



k — kj 




1 — erf 



\V2k 



a 



If one plots this one sees that this instability (c<j^ < 0) corresponds to /c^ < k'j or 
A = 2n/k > \j = 271 /kj. \j is called the Jeans length which sets the scale for 
instabilities. 

The validity of the homogeneity assumption and Jeans effect is legitimate on 
small scales and so stationary stellar systems are generally stable on small scales. 
The clumping of material that begins when A > Aj can sometimes be arrested due 
to some nonlinear effects. 

The Jeans analysis, and significance of the Jeans length, can also be used to 
analyse homogeneous collapsing and expanding systems. The perturbed gravita- 
tional field merely accelerates collapse or decelerates expansion and no Jeans effect 
is invoked. This is the kind of instability that could have led to galaxy formation out 
of the initial homogeneity of the early expanding universe, and so also a compact 
self-gravitating object from a homogeneous soup of bosons. 

1.3 Boson Stars: Nonspherically Symmetric Configura- 
tions 

We have set up boson star configurations in a general 3D code, and used them as a 
test of this code by comparisons to spherically symmetric results obtained using a 
ID code. In the process of stabilizing the code to achieve this a study of coordinate 
conditions was needed. Once stability was achieved the code could be used to study 
the behavior of boson stars under non-spherical perturbations, which we obviously 
could not do in ID. When perturbed in ID a boson star emits scalar radiation and 
loses mass. In 3D we can study another kind of radiation, namely gravitational 
radiation. 
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We now give a brief review on Einstein's equations as a source of gravitational 
radiation and then outline the other works in the field of higher dimensional boson 
star studies. Then we set the tone for our own work in higher dimensions by a 
description of the 3 + 1 ADM split of space-time and introduce the concept of 
extrinsic curvature. A brief description of the tools of numerical relativity and the 
underlying numerical errors that one has to deal with follows. Details of our work, 
and resolution of these problems in the context of our own model, are described in 
the fourth chapter of the thesis. 



1.3.1 Einstein's Equations as a Source of Gravitational Radiation 

Maxwell's equations have radiative solutions that predict the existence of electro- 
magnetic waves. Likewise Einstein's equations have such solutions, giving rise to 
gravitational waves. The construction of laser interferometric detectors like LIGO 
makes the detection of this kind of radiation a real possibility. The theory of grav- 
itational waves is complicated due to the nonlinearity of Einstein's equations. In 
order to simplify the system, and yet get insight into the nature of gravitational 
radiation consider the weak field limit, with the metric just slightly perturbed from 
the Minkowski metric [^, g^^, = rjn^, + h^iy. Since the derivative terms only come 



from the /i^i, the Ricci tensor is 



R 



flU 



Xli,u 



fiu,X 



Oih'), 
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(90) 
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flU 



&" 



dx^ dx" 



h'^ a — h ^^y\ — h i,^^x + \Z\h^lu 



(91) 



From the trace of Einstein's equations R = —SttGT^x. The Einstein equations to 
order h are then 



where 
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^x^^ dx 



-h\ - h\,,x - h\^^x + nV = 16 TT G S-^ 



fiuj 



i^ flU -'' IJLV „ ^/J2^ -^ A- 
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(92) 



(93) 
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Exploiting the gauge invariance of the theory to choose a convenient one 

g'^'^ r ^,y = 0, (harmonic coordinate system) (94) 

gives to order h: 

= \- Iv'" hs,,. + v^" hs.,, - v^-" h,,^s] (95) 



„A(5 



2fcV-^n 



This means -^h^v = \-^h^ pl making (92) a wave equation of the form 

UKu = IQttGS^,. (96) 

In electromagnetism, a locahzed source has no charge flowing in or out of it due 
to charge conservation. The monopole part of the potential of a localized source is 
static, and fields with a harmonic time dependence e*'^* have no monopole terms. 
Hence, electromagnetic radiation is dipolar. Similarly, energy conservation ensures 
the absence of monopolar gravitational radiation. In addition, the power output 
of dipole radiation is related to the second derivative of dipole moment, which 
for gravitational radiation is zero because the first derivative of a mass dipole is 
momentum and the law of conservation of momentum makes its derivative zero. 
That is to say, gravitational radiation is quadrupolar in nature. In spherically 
symmetric spacetimes with spherically symmetric sources, there are no gravitational 
waves. 

This can explicitly be shown. Consider the source free or homogeneous wave 
equation for gravitational waves and write out its plane wave solution 

h^,u{.x) = e^j, exp{ikxx^) + c.c. (97) 

The symmetric polarization tensor should have 3 x 2 + 4 = 10 independent com- 
ponents. However, the gauge invariance introduces 4 conditions on it giving six 
independent components. By a suitable coordinate transformation one can actu- 
ally show that only two are physically meaningful and the remaining components 
vanish. By subjecting the coordinate system to a rotation, one can show that the 
physically meaningful components have helicity ±2 while the others which can be 
set to zero have helicities zero and one. This is analogous to electromagnetism where 
Maxwell's equations give four components of the polarization tensor and using the 
Lorentz gauge this is reduced to three independent components. Without leaving 
the Lorentz gauge a transformation of the vector potential components in terms of 
the derivative of a scalar field allows one to reduce the polarization vector compo- 
nents to two physical components. A rotational transformation shows these have 
helicity ±1. See [|1^ for details. 
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Gravitational Waveform Extraction 

In order to find out how much energy is carried by gravitational radiation, and what 
the nature of this radiation is, one needs to extract the gravitational waveform. In 
addition, one needs different methods in order to compare the energies from the 
different measurements to check that the code is giving reasonable results. We list 
below the various functions that portray these waves in our 3D numerical code: 

• Zerilli Function 

The first step to waveform extraction using Zerilli functions is to write the to- 
tal metric as a sum of a spherically symmetric background and a perturbed metric 
h^y. Following Regge- Wheeler the perturbed metric is written in terms of spher- 
ical harmonics by associating the components with scalars {hoo,hro,hrr), vectors 



{hoe 

be found in 



ho(j,', hrOi 



hrs) and tensors. The details for spherical harmonic expansions can 
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In terms of spherical harmonics, (actually m = Q spherical harmonics since all 
the perturbations we consider are axisymmetric) the perturbed metric hap can be 
written as the matrix 



/ 



N^HoYio 
HiYio 
ho Yiofi 


Hi Yio ho Yiofi 
g H2 Yio hi Yiofi 
hiYio,e R^{k + g£, 













\ 





Yio 

i?2 (^K sinO + G cos 6^) Yio sin 6 ) 

(100) 

where the Regge- Wheeler perturbation functions {Hq, Hi, H2, ho, hi, K, G) are all 
functions of the radial and time coordinates only |^ . Any information about the 
gravitational wave content is contained in these perturbed metric functions so that 
these must be determined. Using the orthonormality relations of spherical harmonics 
these functions can be found from the full metric. For example. 
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Ho = j;pl guYio sine de, 



(101) 



since the / = mode corresponds to the spherical metric N = —jJ'o 9tt sin 9 d9 . 
Similarly, the functions Hi and H2 are extracted. We also get 
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h[ 



Yiofigre sine de, 



(102) 



and ho follows analogously. In order to separate G and K consider 



— (/lee sin^6' - h^ 



G (Yio^ee + sin e cos e Y^ce) 



(103) 



This is not composed of orthonormal spherical harmonics. Therefore one must 
change basis so the ee and (fxp part of the metric can then be written in terms of 
orthonormal tensor spherical harmonics. (See |^ for details.) We find 



cos e Yioe + sin e Yi 



. _ 2 TT r^ {gee sin^ 6* - 
^B^i [I (/ + !)(/ + 2) {I - 1)] sin^ e 
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sin e Yio de. 



(104) 



(105) 



Once one has these perturbed metric components one has all the information about 
the gravitational wave output of the system. However these metric components 
are dependent on gauge choices and, in this form, cannot yield the actual wave 
perturbations. Following [^ one needs to construct gauge invariant quantities from 
the perturbed metric components. When the background metric is Schwarzschild, 
the Zerilli function can be constructed, describing the propagation of even parity 
waves. This is why, when we use this method, we place detectors at the exterior 
region of our system. One defines the Zerilli function p3| : 
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(106) 



where 



k^ = K + r{l-'^)G,r-2 
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(107) 
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The quantity ipz satisfies the wave equation 



0, 



(108) 



where the tortoise coordinate n = r + 2M ln( 2^ — 1) has been used. The scattering 
potential is given by 



VAr) 
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(/(/ + l)-2 + ^)l(/(/ + l)-2 + ^ 
/(Z-l)(/ + 2)(/ + l)^ 



(109) 



• Newman-Penrose Spin Coefficients 

Consider a four dimensional Riemannian space on which a tetrad system of null 
vectors /^, m^, ffi^ and ra^ is introduced. The orthogonal null vectors / and n are 
made by adding and subtracting a space-like unit vector from a time-like unit vector. 
The orthogonal complex null vectors m and ffi are complex conjugates of each other. 
Far from the source the Newman-Penrose scalar 



^A 



Rap-yS n"' m^ rC 771^, 



(110) 



represents an outward propagating wave. For details see |2J, |25|, ^ 

• Bel-Robinson vector 

Constructed in a manner analogous to the Poynting vector in electromagnetism, 
the Bel-Robinson spatial vector is 



pi 



E^^e^^'Es'^. 



;iii) 



It is effective in tracking gravitational radiation. Here Eap and Bajs are the "electric" 
and "magnetic" components of the Riemann tensor. The normal time-like vector 
n" (see ADM section) splits the vacuum four- dimensional tensor Ra^^s into 



^af3 



rp Ryaf35 n" 



(112) 



5, 
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The details of this method are in |2^, p7 |. 
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1.3.2 Axisymmetric Calculations 

Static or equilibrium solutions in the case of spherically symmetric boson stars have 
been described in an earlier subsection. Static solutions also exist for axisymmetric 



configurations (no azimuthal dependence) as shown in p8[. Again, the field itself can 
have a time dependence which does not affect the static nature of the configuration. 
The field is of the form cf) = (f){r, 6')e*'^*, with a metric of the form 

ds^ = -e^- de + e2° (dr^ + r^ dd^) + e^^ r^ sin^ d d^\ (113) 

Here v, a and (3 are functions of r and 9. The three relevant Einstein and scalar 
field equations are set up, and this set of equations define an eigenvalue problem. 
The basic partial differential equations now are of the elliptic kind, and boundary 
conditions are more easily enforced by writing them in integral form using a Green's 
function approach. That is, if for an operator O and functions F and H 

OF{r,e)=H{r,e), (114) 

then F = / dV G{r, , 9, r', 9') H{r', 6') where OG = 5{r - r'). The Green's function 
is then expanded in terms of the Legendre polynomials P2n with equatorial sym- 
metry ensuring that only the even ones appear. The particle number (and mass) 
versus radial derivative of the central field shows marked similarity in profile with 
the spherical case. The maximum mass is 1.05 Mpjn as opposed to 0.633 in the 
spherical case. 

Using this Green's function technique, the analysis is extended to the case of 
rotating boson stars in |]29[. In |3y] a perturbation approach, where perturbations 



were axisymmetric over a spherically symmetric configuration, was considered and it 
was shown that no rotation was possible. However the assumption of axisymmetry 



was also extended to the perturbed boson field. In [^ this limitation was removed. 
The scalar field was itself allowed to have an azimuthal dependence without affecting 
the axisymmetry of spacetime. In [^ this idea was again implemented but with- 



out the drawbacks of the previous attempt. Namely, this time several equilibrium 
configurations were found, not just one or two, for different values of m. In addition 
non- Newtonian configurations were considered as opposed to just Newtonian ones in 
the former case. Also a "spikelike" solution in the energy distribution in the former 
is believed to be caused by inappropriate boundary conditions. This was explained 



in [^. The metric is written in the form 

ds^ = -e'" df + e^" (dr^ + r" dO^) + e^^ r^ sin^ 6 {d^ - u dtf . (115) 



Here z/, a, /3 and u are functions of r and 6. As in |^ the scalar field is taken to be 

-i{at-mf:) (116) 
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where single valuedness as (p ^> ^p + 2tt dictates that m is an integer. This time 
the Green's function associated with the scalar field equation has to be expanded in 
terms of the associated Legendre functions due to this azimuthal dependence in the 
boson field. The expansions are done till order 2n = 12 and numerically the system 
is iterated using estimated values that step by step converge to the correct value. 
The angular momentum 

J = f T^y/^drd9d<p = mN, (117) 

is related to the particle number N = j° (j = conserved current) through the 
azimuthal number m. J = for the m = case as it should be. This means there 
is a conservation law in operation with specific angular momentum j = T^/j^ = m 
being constant for these configurations. This corresponds to the j constant law of 
perfect fiuid stars except that the constant in this case is an integer. The maximum 
mass was Mmax = 1-05, 1.31, and > 2.22 for m = 0, 1, and 2 respectively (in units 
of Mpi/fi). The mass-energy density — Tq for the m = 1 and m = 2 case shows a 
clearly non spherical almost toroidal distribution. The main difference for m = 1 
and m = 2 is the nonvanishing of this function near the origin on the symmetry axis 
in the former case although it vanishes in the latter. 



A rotating boson star configuration with large self-coupling is described in [^ 
This is with a view to having a structure emitting enough gravitational radiation 
as to be capable of being detected by future gravitational wave detectors. It comes 
with the added bonus of a simplified set of equations in which field gradients other 
than those with respect to the azimuthal angle ip can be ignored. As a result of the 
large self-coupling assumption, one can write an effective equation of state. In the 
tail region, where gradients cannot be ignored, the field is itself vanishingly small. 
The Einstein equations are set up and a Green's function method is again employed 
to find solutions. On the left hand side, terms for which the fiat space Green's 
function are known are kept, while the right hand side consists of the other terms. 
As a result the function on the left is the integral of the Green's function times 
the terms on the right side. The Green's function is then expanded in Legendre 
polynomials as described previously. The mass profile and toroidal geometry of 
the density are extracted as per the earlier discussion. The structure of the star is 
determined by the azimuthal quantum number m, the eigenvalue a and the value of 
\f\jw? where A is the self-coupling parameter. As the mass of the star increases its 
radius slightly decreases and the strongest gravity region would be a configuration 
near the maximum mass. These stars would give the largest gravity wave signals. 
From the output gravitational waves, the multipole moment structure of the star 
is revealed. If, from the waves, the mass, spin, mass quadrupole moment, and 
spin octopole moment are determined, then the object could be confirmed to be 
a particular configuration (three quantities completely parametrize it) of a boson 
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star. The maximum mass as a function of the star's spin angular momentum is 
plotted. The multipole moments are encoded in the asymptotic form of the metric 
coefficients and can be determined from comparisons with the same order terms in 
the previously described series expansions. Plots are made of the mass quadrupole 
moment and spin octopole moment as functions of mass and spin. If the gravitational 
wave measurements give four moments, the mass M, spin 5*1, quadrupole moment 
M2, and spin octopole moment 5*3, one could determine from these plots the value of 
the self-coupling parameter and identify the boson star. Plots can also be made of 
the energy of a test particle orbitting and entering a boson star for a particular mass 
boson star, as a function of orbital frequency. One can also make such plots against 
other boson parameters, thereby providing the gravitational energy per frequency 
band. This could serve as a map of the interior of a boson star if there was some 
way of calculating the total energy emitted by measurements from a single detector. 
However, this would involve knowing the exact angular pattern of emitted waves. 

1.3.3 Boson stars in 3D: Perturbation Studies 



The nonradial quasinormal modes of a boson star are described in 33 1. The 



even-parity axisymmetric perturbations of the metric are written using the Regge- 



Wheeler gauge |^0[; with the perturbed metric functions and components of the 
complex boson field having frequency a, that is, a time dependence of the form 
g-«o-t_ rpj^g perturbation equations are then written down for the perturbed Einstein 
and scalar field equations. Outside the star, the background spacetime becomes 
almost Schwarzschild and the scalar fields perturbations decouple from those of the 
metric. The metric perturbation equations then reduce to a system of perturba- 
tion equations for a black hole, and one can then calculate the Zerilli equation as 
discussed above. Since the time dependence is of the form e~*'^*, the second time 
derivative in the Zerilli equation (108) gives a —a'^ipz contribution and one gets the 
radial part of this function Z in the far zone to be of the form 

Z = A^^e-^'^'* + Aoute""'* + ^ (^) ' (^^^) 

since the Zerilli potential Vz is of order \ according to (109). The scalar field 
equations for the two fields reduce to Schrodinger type wave equations and and the 
quasinormal modes for the star are determined by numerically solving the perturbed 
equation for the I = 2 01 quadrupole modes. 

The most important features of the quasinormal modes are that the imaginary 
parts of the frequencies are large and the difference in real parts of the frequencies 
between nearest modes is almost constant. The large imaginary parts of the frequen- 
cies means a damping time scale which is very small compared to other relativistic 
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stars. This is of consequence to our own 3D studies described in chapter 4 and so 
we provide an analysis of the reason for these differences. 

The reason for the differences from other stars conies from the fact that the scalar 
field extends to infinity, as opposed to ordinary stars which have definite surfaces. 
The fact that perfect fluid stars have at least two families of QNM is explained by 
a "two string model" ^A . The star and the spacetime are described by a finite and 
a semi-infinite string respectively. 

For small coupling constants of the connecting string, there are two kinds of 
normal modes. One has a small imaginary part to its eigenfrequency, and for this 
mode the amplitude of the finite string is very large compared to the semi-infinite 
string. This mode is weakly damped. The other mode is strongly damped and has 
an amplitude that is large compared to the finite string. 

Here on the other hand, we have two semi-infinite strings as the scalar field also 
extends to infinity (refer to Fig. 1). 

The displacement in each segment can be expressed as 



y = A exp{iuj{t + x/c)) + B exp(zci;(t — x/c)). 



(119) 



The complex amplitudes A and B characterize the solutions in each segment. At 
X = strings are fastened so that 



vM = ^2(0) = 0. 

Continuity at the position of the spring gives 

yEp{i) = ypx{i), 

ypqil) = VQYil)- 



(120) 



(121) 



There are no incoming waves at oo so Apx = Aqy = 0. In terms of the spring 
constant k the tension T is given by 



T 
T 
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dx ' 
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-kiyp-yg), 
kiyp-yq), 



(122) 



where the ys are displacements of points and segments of the string. Using the 
boundary conditions, and the fact that A and B are constant in each segment yields 
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(123) 
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Figure 1 A schematic diagram of the two string model for a boson star is shown. The 
scalar field, extending to spatial oo, is represented by a semi-infinite string. It is coupled 
to the spacetime (also a semi-infinite string) through a spring of coupling constant k. As 
a result there is no place to store the energy (as in a finite string) and the modes of the 
system are strongly damped. 
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where z = iu l/c. 

If you divide the two equations of (122) by each other you get ApQ = —Aep and 
putting this into either of the equations gives 

ze' = K{e-' -e'), (124) 

where i^ = ^ is the couphng constant. So for small coupling constants, and small 
imaginary parts for the frequencies, we have ze^ = 0. This has the trivial solution 
z = 0, meaning there are no normal modes for the weak damping case. On the other 
hand, for perfect fluid stars the corresponding expression is [^ 



z{e' + e'') = K {e-' - e') ^^—- '-, (125) 

and there are non trivial solutions in the weak damping case. For the strongly 
damped case (e^ subdominant to e~^) both cases yield 

z^Ke~^'. (126) 

If e~^^ is dominant then z must have a real part which is negative. So Re{z) = 
j^^-2Re{z) cQs{Im{2z)) < or cos(/m(22)) < 0. Hence 

z a + i{2n + l)- + ib (127) 

where a is large and positive while b is very small. Writing z in terms oiuj {z = iul/c) 
the eigenfrequencies are approximately 

c f (2n + 1)tt \ . ac ,_,„„, 

By equating the real and imaginary parts of (126) we see that a = Ke^"" (since h is 
small). The presence of a coupling constant, however small, is necessary to ensure 
that we have a nontrivial solution. The smaller the coupling constant the larger 
the damping part of the frequency (a). The imaginary part of the eigenfrequencies 
are equidistantly spaced. These are like the w modes of a perfect fluid star. In 
the perfect fluid case we have a finite string that has no mechanism to radiate its 
energy except through a spring coupling to a semi-infinite string, however weak the 
coupling may be. On the other hand, the two semi-infinite string case allows each 
string to radiate energy excited by an oscillation to infinity without coupling to the 
other string. Thus, there is no place to store energy like in a finite string. Therefore, 
radiation must be rapid and there are no weakly damped modes. 
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1.3.4 Boson Stars in Full GR: Numerical Relativity 

Eventually what one wants is the solution of the Einstein equations with matter 
terms present. This is the crux of our work. Of course, Einstein's equations are 
far too non-linear and complicated to be solved analytically and so we must rely on 
numerical techniques. 

We look at spacetime from the point of view of a Cauchy problem. That is, 
we view the evolution as being from one spacelike hypersurface to the next. In 
order to construct the gravitational field, one solves the initial-value problem and 
then integrates the dynamical equations along trajectories of a prescribed reference 
system. 

Four of Einstein's equations, Gq^ = SttTq^, do not contain second time deriva- 
tives of the metric. The solutions of these equations constitute the initial-value 
problem. The most natural geometrical variables on the hypersurfaces are the met- 
ric of the hypersurface, which is denoted jij or '^Qij, and the extrinsic curvature 

1.3.5 ADM 

In order to develop the equations that describe the system we must understand the 
spacetime that we are studying. In the first place, we must revise our concept of 
time. In Newtonian theory it is absolutely defined. In special relativity time is 
somewhat ambiguous in that the concept of simultaneity is not universal, but by 
specification of an inertial frame the concept is made precise. On the other hand, in 
GR we have to replace the concept of time of an event by the notion of a spacelike 
hypersurface. The entire spacetime is divided into these spacelike hypersurfaces, or 
rather sliced into them, and the parameter separating one hypersurface from the 
other is called "time" . At each event on a given hypersurface, a local Lorentz frame 
exists whose surface of simultaneity coincides locally with the hypersurface. This 
Lorentz frame is one with its 4-velocity orthogonal to the hypersurface. 

Thus to explore spacetime conveniently it is divided into a series of spacelike 
hypersurfaces that are parametrized by successive values of a time parameter t. The 
3-geometry on two faces of a spacetime sandwich are connected by a 4-geometry in 
between that must extremize the action described in |1 . 1| . 

This 3-1-1 split of spacetime, with a hypersurface parametrized by time coordi- 
nate t = constant followed by one parametrized hj t + dt = constant, is convenient 
but this parameter t is not the "proper time" . The 3-geometry of the lower hyper- 
surface is given by 

'-yij(t,x,y, z) dx^ dx\ (129) 

and that of the upper one by 
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•jijit + dt,x,y,z) dx"^ dx\ (130) 

with the lapse of proper time between the two being related through a lapse function 
A^. We have 

lapse of proper time = N(t, x, y, z) dt. (131) 

The spatial coordinates are shifted by a shift vector from one hypersurface to another 

xlgh = <.. + N\t,x,y,z)dt. (132) 

This is part of the gauge or coordinate feedom we have in the theory. 

The 4-geometry connects a point (t, x*) to [t + dt, x^ + dx'^) with proper interval 

ds'^ = 7ij. {dx' + N' dt) {dx^ + N^ dt) - {N dtf = ^g^p dx'' dx^ . (133) 

This gives the 4-metric in terms of the lapse, the shift, and the 3-metric. 

A nongeometrical way of looking at this is to introduce the notion of time through 
a coordinate t = t{x^), /i = 0, 1, 2, 3 and a time flow vector f^. The time flow vector 
is normalized 

t^'\/^t = l. (134) 

Since this vector of time flow is in general not going to be normal to the spacelike 
hypersurfaces, the normal is given in terms of a lapse function 

n^ = -NS/^t. (135) 

Thus, 

t^n^ = -Ar. (136) 

This normal to the spacelike hypersurfaces must be timelike and it is normalized to 
unity so that 

n''n^ = -l. (137) 

Hence 

N = -t^ n. = (n" VJ)'^ = \^ ^ =. (138) 

We now separate the time vector into a spacelike and a timelike part as follows 

tt' = P^^^f -n^n^f = [(5^ + n^n^ - n'^n^)] t'^ = A^^ + A^n^. (139) 
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Here the first term, {5'^ + n^ny)t^ = N^, defines the shift vector. It is perpendicular 
to the timehke component Nn'^, as can be seen by the fact that we have subtracted 
off the timehke component from the whole term to get the first term, and also by 
simply contracting to get 

n^m = {5^,n^t'' + n^n^n,r) = n^t^^^l + (-1)) = 0. (140) 

Since all physics must be independent of coordinates, a system x^ = (a;°,x*) is 
chosen with the spacelike hypersurfaces described by x* and the time hj t = x^. 
The components t^ are then 



The lapse function is 



and the shift 



Also, nfj,n^ = — 1 gives 



Similarly, 



implies A^" = and 



implies 



Thus, 



t^ = (1,0,0,0). (141) 



N = -t'' n^ = -no, (142) 



ni = -NVit = 0. (143) 



n' = ^. (144) 



n^Ar^ = 0, (145) 



N^n>' = 0, (146) 



^ + N,n' = 0. (147) 



n' = — and No = Nib\ where a' = -V (148) 

Choosing 6* = A^* gives 
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nf" = {N-\-N-^N'), (149) 

and 

N^ = {N,N\N,). (150) 

Finally, we can extract the metric components as follows. From 



gij,un^ = riy (151) 

=> goo n° + 9io n' = Uq 
^goo-9^oN' = -N\ 



and 



g>.uN^ = N^ (152) 

^gooN'' + g,oN' = No = N,N' 

one gets 

ds^ = g^^ dx^" dx" = - (A^2 _ j^i jy.) ^^2 ^2Ni dx' dt + -fij dx' dx^ , (153) 

where 7ij is the three metric on the hypersurface. In matrix form 

/ _(iv2_iv*iv,) N, \ 
^-=1 iV. 7J' ^'''^ 

and from gfj,ug'^'^ = S^ one gets the inverse matrix 

• Intrinsic and Extrinsic curvature and Connection Coefficients 

Consider a vector A lying in the hypersurface 

A = eiA\ (156) 

If we parallel transport this vector along a route in the hypersurface and compare 
A at the transported point to the parallel transported A 

V,A = V,(e,- A^) = ej ^ + 'r^, e^ A^ . (157) 
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Here the connection coefficient F is the usual measure of the variation in the basis 
e itself. We started out on the hypersurface but now have a component outside the 
hypersurface A-^^F^eo-n. If one projects ^VA orthogonally onto the hypersurface 
so as to get rid of the component outside the hypersurface then we are only dealing 
with the 3-geometry intrinsic to the hypersurface. Writing the 3-derivative, of some 
vector A that lies on a spacelike hypersurface, intrinsic to the surface and then 
taking the projection in a particular diection gives for covariant components 

8 A • 
ej-3V,A = ^-A™F^,,. (158) 

Here the connection coefficients are those in three dimensions 

^ mji '- mji ^m' VjGj. yio\y) 

The concept of extrinsic curvature deals with the embedding of the 3-geometry 
in an enveloping spacetime. Consider the "normal" vector that stands at a point 
on the hypersurface {V + dV), and compare it to the one parallel tansported to this 
point from a neighboring point V. This can be regarded as the limiting concept of 
the 1-form dn. Hence the 1-form surface lies perpendicular to the vector, and so 
lies along the hypersurface. 

dn = -K(dP). (160) 

The components of the extrinsic curvature Kf are obtained by displacements along 
a particular coordinate direction Cj, so that 

(161) 



(162) 



The change of a vector e^ can be written in terms of its component along the 
normal and along the basis vectors ei i = 1,2,3 and so 

^V^e, = (^V^e,) -1111+ (^V^e,) ■ ek ek = K,, -^ + ^^F^, e^. (163) 





^V,n = -K{e,) = -K^ey 


Then 




Kik 


= Ki gjk = K^ (e, ■ e^) = -(^V,n) ■ e, = -^V,(n ■ e,,) + n ■ ^Viie^ 




= + n ■ 4V.(efe) = (n • e^) ^F^,, = (n ■ eo) ^F\, = n ■ ^V^ei = Kk^. 



Similarly, defining Al = e,- • ^VjA, 



^V.A = A\, e, + K,, A^j-^. (164) 

' n • n) 
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The extrinsic curvature components in terms of the connection coefficients and 
the ADM metric then are 



K^, 



Thus 



'^i;j ^fJ, -'- ij 



K 



1 



'^ 2N 
1 



2N 
1 

2iV 



-no 'T% = -N r°^. = -N {V' To., + V T,,,)) . 



9oi,j + 9oj,i ~ 9ij,o ~ 2 Tkij N 



dx^ dx^ dt 
dgij 



2Tk.,N'' 



iV,|, + iV,|, 



dt 



(165) 



(166) 



• The 3 + 1 Gauge Freedom 

One sees, therefore, that the spacetime metric g naturally separates into the 
six jij (symmetric), the three shift terms iV*, and the lapse N relative to a given 
foliation and choice of time vector. The only second-time-derivative terms in the 
Einstein equations are those of the ji/s and these are the dynamical degrees of 
freedom. However the initial-value problem, that we have hitherto not discussed 
and which we shall address in chapter 4, shows that the 7ij need only be known up 
to an overall conformal factor 0"^, with determined by constraints. Since 



7ij 



^Y' %• 




(167) 



(7)^/^^,, 7 



one can regard the separation of the configuration coordinates 'jij into the square 
root of its determinant ^/^ and the conformal metric 7jj. Since tr K turns out to be 
a canonically conjugate variable to ^7 one can as well regard the six {'yij,trK) as 
configuration variables. However, we have four constraint equations and hence only 
two dynamical degrees of freedom. Assuming the constraints have been satisfied we 
then have four gauge degrees of freedom that we impose as four conditions on the 
velocities dtitrK) and dt'^ij. These lead to equations on the lapse and the shift [ ^ . 
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• The 3 + 1 Equations 

Consider the projection n^n^T^y = n^riuT'^'^, where n is the normal vector de- 
fined in the earher subsection. This is n^riQ Tq g^^ = n^riQ Tq which is the energy 
density p = —Tq. The Hamiltonian constraint equation is determined from 



n^n^G^^" 



til p. 
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The projection —n^P^iTi^fj, where (139) defines the projection operator Pi, is 
—vy5fTyfj — rfnin"Ty(j. Since n^ = it simphfies to —rfTyi = —vyTi^go^ = — noTj°, 
which is just the momentum density J'i. This gives three momentum constraint 
equations 



-n-^p.^a 



8 71 J,. 



(169) 



Finally, the projection Pi^Pj^T^^ = Tij = Sij provides six evolution equations. The 
extrinsic curvature relation (166) in the 3 + 1 approach, guides the evolution of the 
metric through six more equations. 

We can write G^^ in terms of the Ricci tensor i?^j, and Ricci scalar R, which we 
can write in terms of the Fs' themselves. These are then written in terms of the 
extrinsic curvature, to get the constraint and evolution equations: 



^R - K] Kl + K' 



,71 p, 



Hamiltonian Constraint 



(170) 



K,'\,-K\ 



871 Ji, 



Momentum Constraint 



(171) 



and evolution equations 



1 

'n 



d_ 



K/ - N\k Kj^ + N^^j Kk' + N^' 



Stt 



Si -^iS-p) 6/ 



+ -'R/ + KK/ 



(172) 



and 



0, "jij 

IT 



+ N\nkJ + N\,^k^ = -2NK 



iji 



(173) 



where S = Si and K = Kl. Instead of six second-order evolution equations for 7jj, 
by defining extrinsic curvature in terms of the first time derivative of the metric we 
have 12 first-order equations. By letting j = i and summing over all i in (172) we 
get an evolution equation for tr K 
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^lL^ = ^N\'.i + N \r + K^-Att(3p-S)\. (174) 

at L -I 

Data on the initial hypersurface satisfying the 4 constraint equations evolves 
to the next hypersurface, according to the 12 evolution equations. On the next 
hypersurface, but for the shortcomings of numerical approximation schemes, the 
data should consistently satisfy the constraint equations. 

• Coordinate Conditions 

Part of the gauge freedom we have is that we can impose four conditions on 
dttr K and dttr^ij as discussed before. This gives us equations for the lapse and 
the shift. We will discuss the various coordinate conditions we use on our lapse and 
shift, and the specific advantages of these in chapters 2, 3, and 4 as is pertinent to 
the problem being considered. 

• Finite Differencing Errors 

There are three kinds of partial differential equations that one has to deal with. 
Our evolution equations are typically hyperbolic equations (or wave equations) and 
the constraint equations are elliptic equations. In addition, there are parabolic 
equations, and parabolic terms may be added to our equations as part of developing 
a numerical algorithm. Elliptic equations represent equilibrium and the constraint 
equations are of this form. In addition, when we use maximal slicing (described 
in chapter 4) we have an elliptic equation for the lapse that has to be solved on 
each time slice. This is in fact the most time consuming part of the numerics. The 
constraint equations are solved on the initial timestep and then monitored as the 
code evolves, on each time step, as indicators of the accuracy of the code. We 
circumvent the 3D constraint problem, for spherical configurations, by getting data 
from the ID code and putting it on a 3D grid. We are dealing with highly nonlinear 
equations but to get an idea of the nature of these equations, we write them down 
in the linear case. 

dtdtcp — (? dxdx4> = hyperbolic (wave) equation, (175) 

dydyCp + dxdxCJ) = elliptic equation, 
dtcf) — Ddxdx(p = parabolic (diffusion) equation. 

The factor D is a diffusion coefficient. 

We have already shown how the linearized Einstein's equations give rise to the 
gravitational wave equations. In fact the full evolution equations are like the wave 



equation I^Tf. By defining the extrinsic curvature in terms of the time derivative 
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of the metric we evolve a set of first order equations. So a prototypical hyperbolic 
equation to consider is of the form 

(9t0 + t; (9^0 = 0. (176) 

We use finite differencing techniques to solve these equations. 

The "spacetime" for the numerical algorithm consists of a discrete set of points 
and not a continuum; the grid spacing, and timesteps, being dictated by the re- 
sources of memory and real time available to us. The derivatives in the PDE are 
replaced by finite difference approximations. In doing so one introduces spurious 
solutions of the finite difference equations that have no bearing to the solutions of 
the actual equation we are solving. There could be several different ways to finite 
difference the problem, each with its own solution. Choosing a finite differencing 
scheme is, therefore, of great importance. 

One obvious kind of error is the Truncation Error. We typically use second order 
accurate schemes for finite differencing in the 3D code and up to 6th order in the 
ID code. There can also be Propagation Errors. For example if one writes the 
prototypical hyperbolic equation as 



iji+i 



G0", (177) 



where G is an operator connecting the function on time slice n with that on time 
slice n + 1, one can have propagation errors. If the solution of the finite differenced 
equation exhibits a wave behavior that is not shared by the PDE it represents, the 
solution can grow without bound ii \G\ > 1 and it will dissipate away if |G| < 1. 
If G = 1 then the wave amplitude is propagated without change as the physical 
solution. For example, a forward-time centered-space derivative scheme can become 
unstable. We discuss this in chapter 3 under our leapfrog method discussion for the 
numerical code. In the staggered leapfrog method, one uses second-order accuracy 
in space and time and |G| = 1. 

In nonlinear hyperbolic equations, we often see grid point to grid point instabil- 
ities due to our using a leapfrog scheme. In a leapfrog scheme, one uses centered 
time and space derivatives. As a result the time levels in the time derivative terms 
"leapfrog" over time levels in space derivative terms. For example, in our prototyp- 
ical hyperbolic equation 

ur'-ur' = ^{u]^,-nU)- (178) 

This kind of method for a non-linear equation can result in odd and even mesh 
points totally decoupling, and results in mesh drifting with odd and even grid points 
evolving in their own separate ways. For large gradients, this can result in grid point 
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to grid point instabilities that can invalidate the code. To prevent this, one needs 
to connect odd and even grid points through a viscosity term. This is a term of the 
form D ("Uj+i — 2-u" + Mj.i) where D is a small parameter. We have had to use a 
very tiny term of this form in our 3-D code. Putting this kind of term everywhere 
usually invalidates our code at our boundaries so that we have to exponentially 
damp it out in the exterior of the star. The steep gradients occur only near the 
center of the star, so we never need diffusion anywhere else anyway. Typically we 
add diffusion in the central l/3rd of the grid and then damp it out. 
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Chapter 2 

Equilibrium Bosonic Objects in Spherical 
Symmetry 

We now go on to continue the analysis of the Einstein equations with matter sources 
in a spherically symmetric environment. The general spherically symmetric metric 
is of the form 

ds'^ = -{N^ - N'' Nr) df + 2Nrdtdr + g'^ dr^ + r^ {dO^ + sin^ 6 d(j)^), (179) 

where the metric functions A^^ and g^, and the shift function Nr, are only functions 
of r and t. The angular metric functions are gge = r^ and f?^^ = r^ sin^ 9. 

We study the stability and dynamics of such a system using a numerical code. 
The success of a numerical code relies on the choice of coordinates. A general code 
with general metric functions and x, y, z coordinates would have to perform far too 
many algebraic operations to calculate all the Ricci components. A good coordinate 
system takes advantage of the symmetries of the system and avoids physical singu- 
larities by slowing down evolution in the spatial region near the singularity. Also, 
it should itself not develop any coordinate singularities. 

The lapse function determines the proper time separation between nearby space- 
like hypersurfaces and determines the time coordinate. The choice of lapse is often 
times determined by a condition on the extrinsic curvature (refer to |1.3.5| ). Our ID 
code uses a polar slicing condition Kg + Kt = [Q . In a spherically symmetric 



spacetime with no A^^ and A^"^ terms this translates to 



'-N 



reeg^' + T'-^^g'^^ N, (180) 
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Thus the shift Nr is zero. 

The polar slicing condition in a spherically symmetric spacetime is very singular- 
ity avoiding, matching up with Schwarzschild coordinates. In the spherical collapse 
problem (in this case when a boson star is collapsing to a black hole), the evolu- 
tion slows down so that not even an apparent horizon forms as the lapse collapses 
(A^ ^ 1 — 2M/2M = 0). Thus when an apparent horizon is approached, the nu- 
merical code breaks down as the lapse collapses and the radial metric rises sharply. 
This is because sharpening gradients invalidate the finite difference approximations. 
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We now outline the features of the various studies we have conducted in spherical 
symmetry. 

We first calculate the connection coefficients from (4). 

(181) 
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Here prime refers to a derivative with respect to r and rfot to a derivative with 
respect to t. 

The rest of the coefficients are zero (other than the symmetric exchange of the 
lower two indices) as there is no variation with 6 or (p. These are the connection 
coefficients we must now use to calculate the Ricci tensor and the Ricci scalar, R^^, 
and R. 

ID pA pA pA I pA pp pA pp (^flO\ 

From the Riemann tensor we calculate the Ricci tensor and Ricci scalar to give 



us 
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G^, = R,,--g>'''R. (183) 
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This gives 
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Since Goo = SttToo according to Einstein's equations, we get an equation for g' 
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The rr component of the Einstein equation is Grr = SirTrr- This gives us an equation 
for A^' which is 



N' = -N 
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Similarly from the tr component of Einstein's equations we get 

g = A^TGrgTtr■ 



'l%9) 



These are the equations in GR. Equation (187) is the Hamiltonian constraint equa- 
tion and is used to monitor the accuracy of the code. Equation (188) is used to 
calculate the lapse on each time step and equation (189) is an evolution equation 
for the radial metric which we use to get the radial metric on a given time step from 
its value on a previous time step. In Brans-Dicke theory, which we will discuss in 
the next section, the form of the field equations changes and the StiT^u terms in the 
above equations would have to be replaced with those plus additional terms. 

2.1 Alternative Theories of Gravity 

We have also studied boson stars in Brans-Dicke theory, and now give an overview of 
the history of these theories. We then derive the equations for spherically symmetric 
boson stars in a Brans-Dicke theory. By doing this we are also able to get the 
equations in GR by setting the BD parameter cj = oo in these equations. 

Alternative theories of gravity, like Brans-Dicke and Scalar-Tensor theories, were 
first propounded to explain the coincidence in values of large numbers made from 
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dimensionless combinations of physical and cosmological constants. The idea as 



originated by Dirac and outhned by Jordan ||39[ was as follows: 

There are six constants that determine the large scale structure of the universe 
namely 

• 1) The velocity of light c. 

• 2) The Gravitational constant k = ^^. 

• 3) The age of the universe A. 

• 4) The Hubble effect red shift to order 1/r, is, Ho = ^^. Here the observed 
nebula is at distance r from the observer and its observed red shift is AA. 

• 5) The mass density throughout the universe fi. 

• 6) the radius of the universe R. 

Amazingly all the dimensionless quantities HqA, ^, and kjjic^A^ are of order 
unity. The first of these relations on substituting for Hq shows that when the age 
of the universe was small the expansion must have been rapid. The second relation 
suggests the radius of the universe is increasing proportional to its age. The last 
relation in the context of a homogeneous and isotropic closed universe can be written 
as 

kfic^A^r^k — R^ = ——^l. 190 

it'' R 

Since the radius of the universe is not fixed it means that either k, M or both are 
varying. 

Turning now to constants derived from microphysics: consider the ratio of the 
radius R of the universe to the elementary length / (the electronic radius e'^/rriQC^). 
This is a number of the same order of magnitude (10^*^) as the ratio of the electro- 
magnetic to the gravitational force, and since R/l increases proportional to the age 
of the universe, it suggests that the gravitational constant might not be constant 
but might be inversely proportional to the age of the universe. 

Assuming one accepts that dimensionless quantities made from physical and 
cosmological constants are of the same order, then the gravitational constant is cos- 
mo logically not after all a constant. Under this premise a modified relativistic theory 
of gravitation compatible with Mach's principle was developed |^^ . In Newtonian 
and earlier theories space itself is portrayed as having its own physical structure and 
properties. Mach's view on the other hand was that the physical properties of space 
were intimately connected with the matter therein. The only meaningful motion 
of a particle would then be motion relative to another. GR partially accomodates 
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to this point of view, with spatial geometries being affected by mass distributions. 
However they are not uniquely specified by those distributions. In order to fully 
incorporate Mach's principle, boundary conditions would have to be imposed on the 
field equations, eliminating solutions with no mass present. 

As seen earlier, the gravitational constant, under the premises of the large num- 
bers hypothesis, is dependent on mass distribution in a uniformly expanding universe 
through 

GM , , 

^~l. (191) 

Here M is the finite mass of the visible universe, assumed to be of finite radius 
R. This relation suggests that either the M to R ratio is constant or the locally 
observed gravitational radius (G/c^) is determined by the mass distribution about 
the point in question. The former might come out of some boundary conditions 
on the GR field equations. The latter however is not compatible with the "strong 
equivalence principle" and GR. 

As per Mach's principle, the motion of a test particle is only meaningful in rela- 
tion to the rest of the matter in the universe, and a laboratory accelerated relative 
to distant matter must experience inertial forces equivalent to a laboratory that is 
fixed due to the existence of distant accelerated matter. If the inertial reactions are 
interpreted as gravitational forces due to distant accelerated matter then locally ob- 
served values of a particle's inertial mass could depend upon the matter distribution 
around it. However there is no way to compare the mass of a particle at one point 
with the mass of another particle at a different point. Mass ratios can be compared 
at two different points but not masses. Using the gravitational constant we get a 
characteristic mass the Planck mass. 




= 2.16 X 10"^^, (192) 

with 

me f ^y ^ 5 X 10-23, (193) 

providing a way to compare an electron's mass at different spacetime points. Like- 
wise the notion of h and c being the same at all spacetime points is meaningless. 
However h and c can be defined to be constant without ambiguity if m or G or both 
are allowed to vary from position to position. (Earlier we discussed the possibility of 
variation in time.) The formal structure of the theory would be different depending 
on whether the mass changes or not. 
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In GR the representation uses units so that atoms have physical properties inde- 
pendent of location. Assuming such a choice is possible the theory is simpler if the 
gravitational constant is the one allowed to vary with inertial masses of elementary 
particles being constant. In GR, the strong equivalence principle states that the 
laws of physics, and their numerical content (like dimensionless physical constants) 
are independent of the location of the laboratory in spacetime. Mach's principle 
as outlined above is only compatible with the "weak equivalence principle". This 
regards local equivalence of gravitational forces and accelerations and is the only 
one experimentally supported by the Eotvos experiment. 

The Brans-Dicke theory has its foundation on the above arguments. A gener- 
alization of GR is made in which the variation of the gravitational "constant" is 
effected by writing it as a function of a scalar field which is coupled to the mass 
density of the universe. The contracted metric tensor is a constant and is uninter- 
esting. The scalar curvature and other scalars formed from the curvature tensor 
contain gradients of the metric tensor components and fall off faster than 1/r from 
a mass source. Hence these are determined by nearby rather than distant matter. 
Some new scalar rather than those provided by GR must therefore be introduced 
which would determine the local value of the gravitational field. This scalar field 
satisfies a scalar field equation of the form 

nv = ^^^T%j^ (194) 

where A is the coupling constant tying the field to matter sources. The average value 
of the field can be roughly estimated by taking the universe to be a uniform sphere 
of density p ~ 10~'^^g/cm^ and radius R ~ lO'^^cm. This gives an average value 
of the field (~ \pR^) of roughly 1/G (units c = 1) if A ~ 1. These considerations 
suggest an equation of the form 



G --^- 



T^v^ + T^j, . (195) 



However one does not want to give up the successes of the equivalence principle, 
namely the equality of gravitational and inertial masses, so the BD parameter should 
not enter into the equations of motion of the particles leaving T^ ^-i, = and of the 
same form as before. 

One starts from the usual variational principle of GR 



= 6 



R H -, — Lm 



d" 



-g d^x, (196) 



(with scalar curvature R and matter Lagrangian Lm including all nongravitational 
fields). Dividing through by G, and replacing the inverse of G by a scalar field (/? 
gives a scalar field equation: 
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= 6 



ipR^ 7-Lm - ^ {^,f, V^'^/v^) 



c4 



[-g)-^ d^x. (197) 



Note that a Lagrangian density for the scalar field has also been added. In the 
original Brans-Dicke theory uj is constant. Writing ip'^ = g^^(p^v and finding the 
variations as we did in the Introduction we get, using the traditional c = 1 units 

-R/ii. - ^fi'M'^-R = 8 7rv9"-^T^!. + — Up^^Lp^i,- -g^yLp^^Lp''^] (198) 

Here T^iy is the usual source term from (29). All the right hand terms have been 
multiplied by the factor of G or p~^ in this case. 

2.1.1 Boson Stars in Brans— Dicke Theory: Spherical Symmetry 

We now turn to the description of Boson Stars in Brans-Dicke theory and derive the 
equations that describe the system. We start out with this model because Boson 
Stars in GR can be derived as a special case of these and the equations for the case 
of massless scalar field halos follow as a further specialization. These are the three 
spherically symmetric models that we have studied using our ID code. In fact we 
start with a general scalar tensor case where the parameter u is not constant but 
may in fact be a function of the scalar and/or tensor fields. 

The matter Lagrangian in this case is one from which the relativistic scalar field 
or Klein-Gordon equation is derived and is of the form 

LM = -i^^^9^0*9.0-im2|0|2-iA|0|^ (199) 

where is the boson field and the A |0| term is a self-coupling term. This is 
the simplest way to include a self-interaction while ensuring stability as well as 
renormalizability. 

The action for the system is given by 



Sj = -^ J d^x ^f^j UpjR ^d^pjd^pjj + J d^x ^[^jLm, 



(200) 



where the subscript J refers to the fact that we are in the physical or Jordan frame. 
By a "physical frame" we mean a frame in which all fields give rise to positive- 
definite energy density. In the Jordan Frame gravity is entirely described by the 
metric g^^. There is also another frame, related to this one by a conformal transfor- 
mation, the Einstein Frame. In this frame the scalar tensor field acts like an external 
matter field which is the source for the metric. Hence, in the Einstein Frame, where 
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the equations are like those in GR, gravity is not represented by the metric alone 
but by the scalar-tensor field as well. Thus scalar-tensor theories and GR may be 
mathematically similar, but are physically dissimilar. The field variables Qfj^iy and 
(fj define the Jordan Frame. By assumption the matter couples to the metric and 
does not couple to the gravitational scalar field. Thus the Jordan Frame is assumed 



to be the physical frame |41 



Varying the action with respect to the metric gives, as seen in the previous 
subsection, the gravitational equation 

G,. = —T,^j + —^ (d,^j d.^j - ^ d,^j 9Vj) (201) 

ipj if/ \ 1 J 



+ — (</';/.;i.j-^Mi'j</';A'^j) 



Varying the action with respect to the boson field, gives us, according to the Euler- 
Lagrange equations 

5a^ = ^' (202) 

d{dx<p) d(f) 

the scalar field equation 

g'^'^j dxd,<P 6\-m^<p-X |0p 0* = 0, (203) 

or 

0;A'^-mV-A|0p0* = 0. (204) 

Varying the action with respect to the ST field (f gives 

—Vj'\,x - 4 d'vj dxipj + Rj + -^^ d,vj d^^j = 0. (205) 

Taking the trace of (201) gives 

i? = -— T + 4a.^aV + ^. (206) 

Also, from the form of the matter Lagrangian, we have that the stress energy tensor 
in this frame is 

T,. = I {d,<P* d,<p + c.c) - \ g^\j [g^\j d,<P* 9.0 + m^ 10^ + i A 10^) . (207) 
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However for convenience we make a conformal transformation to the Einstein 
Frame keeping in mind that when we solve the equations we must revert back to 
the Jordan Frame when interpreting the results physically. 

The transformation from one frame to the other is given by 



We define a constant G* through 
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(208) 
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(211) 



Here the quantities without subscript J refer to quantities in the Einstein frame. 
Going back to the action in the Jordan frame (200) we get from 



Rj = e-''^ [R - Gg^-'V^V.a - 6 g^" (V^a) (V.a)] , 
and the other transformations 
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or 






(215) 
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Since V^V^a is a total divergence the integral vanishes and we get the action in the 
Einstein frame 

Se = -^ [ d^x y^ (R - 2 g''" V^if V,v^) (216) 
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2 '"' 4 
Again using the variational principle described in the earlier subsection we get 

G^. = 8 TT a T^. + 9^9> - dr^ 9> (7^„ (217) 

We get the stress energy tensor from the matter Lagrangian term by variation with 
respect to the metric as seen in the previous subsection so that 



T, 
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-^e^"^^) (9,0* 9"0 + e^"^^) 



(218) 



m^ 0* + 2 \/(0*0) 



9n.vi 



where the first set of terms comes from variations of Lm and the second set from 
variations of \J —g with respect to the metric. Variation with respect to the Brans- 
Dicke field gives from the E-L equation 



-2 g^'' Va^; V,v^ - 2 g^-Vx St V^y^ 



(219) 



-16 7re^'^2a^'^'^ [ - d^c/) d^(l)* ] + 167re^Ma -m"|0r--A 



So that 

VaV^V? = -Ana 
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(220) 



:22i) 



Va V^V^ = -4 TT a T, 

where T is the trace of the stress-energy tensor. For the scalar field the E-L equa- 
tions gives 

V^ (e'" ^^0) = e^" (-m^ - A |0|' 0) , (222) 

giving 

5A,„ „2a f _2 J, \ I ^|2 



V^V^0 + 2 a 9a0 9V e''^ -m' - A 



(223) 
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2.1.2 BD Stress Energy Tensor: Equilibrium Conditions 

From (218) we see that the relevant components of the stress-energy tensor for the 
spherically symmetric problem are given by 



-00 
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^|9.0p + |0p + iV2e^«L^|0p + ^|0r' 
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(224) 
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Note that we had been using units with h = c = 1 so that distance and time 
were in units of l/mass. The boson field has an underlying oscillation frequency uje 
for an equilibrium boson star. If we look at the form of the stress energy tensor we 
notice that every term with the boson field or its derivative is multiplied by 0* or 
its derivative, indicating that a time dependence in the field of the form e*"^-^* does 
not introduce any time dependence in any physical quantity. Hence even though the 
field has this time dependence the system can still be in equilibrium provided the 
metric is not time dependent. This underlying frequency is crucial for a boson star 
configuration to be held together. Gravity, ever attractive, is trying to cause the 
system to collapse and to prevent this collapse the dispersion of the wave equation 
must be used to balance it. For a given configuration a very specific frequency of 
the wave will exactly balance the gravitational effects. We have used the symbol ue 
to distinguish it from the scalar tensor parameter u. 

We now make a change of coordinates to a set of dimensionless ones: 



r = mroid, t = ujeI 



old, 



We define (f){r,t) = <l>(rje 
We now use equations 
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|^184)-(189), with the modification that in the BD case all 



factors of ^nGT^i, be replaced by the right hand side of (201). For the equilibrium 
configurations we take the BD field to have no time dependence. So in dimensionless 
coordinates we have the equations for an equilibrium boson star, namely the scalar 
field equation for the BD field ((220) with ip replaced by (Pbd), plus the scalar field 
equation for the boson field, and the tt and the rr components of the gravitational 
field equation 
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(226) 
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(228) 
(229) 

(230) 



We have just derived the equations of the equihbrium BD boson star system and 
we will now discuss the details of the configurations and code. We note that for 
the case of boson stars in Einstein Gravity we just need to set a = in the above 
equation as there is no BD field involved. Finally, we make the further modification 
m = to get the equations for configurations with a massless boson field. 



2.1.3 Constraints on the value of u; 



The perihelion shift of Mercury provided the first constraint on lu [jT9[. The perihe- 
lion shift is a measure of orbital precession of a planet that does not follow a closed 
ellipse. If it did the change in azimuthal angle at the end of each revolution as it 
went from its perihelion position to aphelion and back, would be 27r. The deviation 
from 27T is a measure of the perihelion shift. 

An isotropic static metric {ds"^ = g^^dx^dx^ time independent and all spatial 
terms made from scalar products of x and dx) can be chosen to be written in the 
standard form ds^ = —B{r)dt'^ + A{r)dr'^ + r'^dQ'^. For a static spherically symmetric 
body like the sun, we expand the metric in factors of MG/r 



, 2 A MG ^M^G^ \ ,2 / MG 
ds^ = -ll + a + f3 ^- + ... rfr + 1 + 7 + ... 



dr^ + r^dn'^. 

(231) 



For the gravitational case in the weak field limit Einstein's equations deliver the 
Schwarzschild metric with a = —2, (3 = and 7 = 2. For the Brans-Dicke case the 
equations give a = —2, /? = 2 — 7 and 7 = 2^^^. Once the form of the metric is 
determined the geodesic equations can be written out for a particle in this metric, 
as, for example a planet in the Sun's metric. Since the field is isotropic we can 
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choose to remove 6 dependences and confine calculations to the ^ = ^ plane. We 
thus arrive at a formula for ^ with two constants representing energy and angular 

momentum. At the major and minor axes of the ellipse, where ^ = 0, one can 
solve the equation to determine the two constants. Integrating from the minor to 
the major axis and multiplying by 2 gives the change in cj) in one revolution and 
subtracting 27r provides the perihelion shift. This effect is maximum for Mercury, 
the closest planet to the Sun. The GR prediction is a perihelion shift of about 43.03 
seconds per century. This is in close agreement with observation. However one of 
the assumptions of the calculation was that the Sun was a perfect sphere, while in 
fact it is slightly oblate. Dicke and Goldenberg |^ photoelectrically scanned the 
solar disk and reported that the Sun's polar diameter was shorter than its equatorial 
diameter by some (0.005 ±0.0007)%. Assuming the assymetry persisted throughout 
the Sun, the imphed quadrupole moment could account for an extra precession of 
3.4 seconds per century. Thus only only 39.6 seconds could be a general relativistic 
effect. Thus Dicke and Goldberg claimed that the GR prediction was really in 8% 
disagreement with observation. For BD theory the value is {"iu + 4)/{'iuj + 6) of the 
GR value so a bound on uj would be given by, 

1^ (108/100) > 1 (232) 

This put a lower limit uj > 6.4. Historically, this was the crucial experiment that 
brought the Brans-Dicke theory into prominence. 

Recent experiments ^, ^ place a much tighter constraint lubd > 500. Note 



the LJ = oo limit of the BD theory is GR. 



2.2 Boson Halos as Dark Matter 



One of the boson models we have studied in ID is made from a massless complex 
scalar field. We have studied the evolution and the formation of such a system [|^ . 
The model had been proposed as a dark halo model [|^, ^ 



As discussed earlier, the flatness of the galactic rotation curves indicates the 
presence of a dark halo. If there was no dark halo, then beyond edge of the luminous 
core, the rotation curves should follow a Keplerian v ~ J^ fall off. Instead they 
are flat, that is w ~ constant which implies M{r) = vf^^^^r. This suggests a linear 
relation between mass functions of galaxies with distance and a density distribution 
p ~ 1/r^. A model for the mass of this functional form, but that does not fall inward 
to the center of the galaxy and form black holes, is needed. 

Attempts to fit the galactic rotation curve using the boson star model have been 
made in ^, ^. While the ground state boson stars have a mass falloff which is 
too fast to explain rotation curves, for more than three nodes the rotation curves 
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resemble observed ones. Since excited states of boson stars are inherently unstable, 
one justifies these models on the grounds that the time scale of instability is large. 
There is a limit on the number of nodes, since the mass of the system increases with 
number of nodes, and Q > 1 ii more than six nodes were present. The best fits are 
then with four or five nodes. The similarity in "quantum numbers" for all galaxies 
is explained on the basis that on an astronomical time scale all galaxies were formed 
simultaneously, the dark matter having collapsed from some more dispersed state of 
higher quantum number to the present one. The absence of lower quantum numbers 
is speculated to be due to higher energy level states being closer together while 
lower energy levels are more separated. Since the energy dissipation mechanism is 
inefficient evolution to lower states would be harder. The calculated density profile 
is p ~ J--1-6 Tj-^ig ig ghghtly increasing compared to the fiat curve prediction of 
p ^ 1/r^. By including a repulsive self-coupling term one can increase the time 
scale of stability. Also one has a wider range of mass and self-coupling to fit curves 
with. However regardless of the value of self-coupling, the rotation curve slightly 
increases compared to the prediction with p ^ r~^'^ . These attempts to fit unstable 
configurations of boson stars to explain galactic rotation curves could be considered 
rather fanciful. A more natural choice seems to involve massless scalar fields. 

The massless scalar field model is governed by equations that can be derived 
from the BD equations above by setting m = and a = 0. To get an idea of the 
nature of the solutions we consider the M — i> limit Newtonian system described 
in 



46| , ^l\ . The scalar field equation is then of the form 



^^2V'-V2V' = 0. (233) 

For a field ip = cr(r)e*'^* this equation has a Bessel solution 

a{r)=A — ^ ^. (234) 

ur 

Regularity at the origin rules out the "cosine" solution. (This provides an ap- 
proximate solution to the low density system. We have studied the evolution of 
configurations of this form in our evolution code to see if they settle down to stable 
configurations. This is described in the next chapter.) 
The density p = — Tq is for this case given by 



A^ 



sin(2co'r) sin^(u;r) 



(235) 



For small r this behaves like a constant and for large r has a 1/r^ dependence making 
it a viable halo model. The rotational velocity profile f|/eu>t = M{r)/r ~ A^. The 
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mass scales as - while the energy density scales as uP' . The central amplitude of 
the scalar field, namely the parameter A of the Newtonian solutions, determines 
the orbital velocity while the scalar field frequency varies the mass and density near 
the galactic center. The best fits to the rotation curves are found if an almost 
dominating scalar field halo is added to low luminosity galaxies. The energy density 
decreases in leading order as J^ jr^iip' and must eventually attain an equilibrium 
with some other matter at the same pressure. Alternatively the radius of these 
entities may be able to estimate the size of the cosmological constant A. 

We have investigated the nature of these solutions in full GR as well as studied 
the formation of these entities. We will describe the equilibrium configurations in 
the next section. The results of our investigation of their stability and formation 
are detailed in chapter three of this thesis. The main question that remains is the 
basis for truncation of this halo at some outermost value of radius. 

2.3 Equilibrium Configurations 

We now describe the results of numerically integrating the equilibrium equations 
for the three spherically symmetric cases. Boson stars in GR and i?D, and Boson 
Halos. We show the mass profiles determined from our code and in the next chapter 
we discuss the results of dynamical evolutions. In all cases the equilibrium profiles 
have static metrics although the fields themselves have an e*'^^* dependence which 
leaves all physical quantities time independent. 

2.3.1 GR Case 

The equations in the dimensionless coordinates defined in (225) for a self-coupled 
equilibrium boson star are given by (227)-(229) with parameter a set to zero (and 
no BD scalar field equation). Note that there is a self-coupling term which we have 
already motivated in the introduction. For very high values of the self-coupling 



parameter, it has been shown in |T3[ that the field varies very slowly over a very 
large radius. As a result of this one can ignore the field derivative terms in the field 
and metric equations as compared to the field itself. In fact in the A — > oo limit one 
can write an equation of state with the field effectively vanishing at the boundary. 
Thus the perfect fluid theorems of stellar stability might be used for very high values 



of A |T^ where an approximate equation of state might be written. These theorems. 



however, cannot apply for smaller values of self-coupling. We choose to study the 
stability of these stars numerically ||3^, |50| . 



The system of equations describing the equilibrium configurations (static metrics 
although the scalar field itself has an e^^^^ dependence), form an eigenvalue problem 
that we have to solve in order to get the equilibrium boson star solution. This is 
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because we have an overdetermined system. We specify the field at the origin. The 
radial metric at the origin is by regularity required to be unity. The field derivative 
and derivatives of the metrics must be zero at the origin by conditions of spherical 
symmetry and absence of divergences at the origin. This gives five conditions with 
four first order equations. As shown before we have the freedom to rescale the lapse 
with the underlying frequency of the boson field. Thus the eigenvalue giving the 
value of the lapse at the origin is physically a way to calculate the frequency of the 
field such that the system has the exact amount of dispersion required to prevent 
gravitational collapse. For a given value of central field density o"(0) the different 
eigenvalues correspond to different numbers of nodes. 

We also have a boundary condition on the boson field to vanish at cx). We 
therefore use a shooting method to determine the eigenvalues. This is done by 
making a guess of the uj rescaled lapse at the origin, and then using this guess to 
integrate the equations. If one is looking for a ground state solution one knows 
the field must have no nodes. It also should be monotonically decreasing and must 
vanish at cxo. So at some stage with this guess the field will either go negative or 
will develop a positive slope. The trick is to bisect between two values one of which 
gives a positive slope and the other of which causes the field to become negative 
and to continue this process as far as the machine permits so that the field falls off 
giving us the solution to greater and greater radii until the machine limit is reached. 
In the one node case a guess value will be used until the field crosses zero a second 
time, or if the slope changes sign a second time and so on. A fourth order Runge 
Kutta technique [^ is used to integrate the equations. 

The mass profile of boson stars is a much discussed feature which we have already 
commented on. However, for the sake of completeness, we show the ground state 
mass profiles as determined from our equilibrium code in Fig. 2. Clearly the mass 
increases with increased self-coupling. 

The mass profiles of excited state boson stars are similar to ground state stars. 
Fig. 3 shows the mass versus central density curves for ground, first, and second 
excited states of boson stars without self coupling. The maximum mass increases 
with the number of nodes as expected. 

• Runge Kutta (RK) Method 

Whenever we have ordinary differential equations of any order we are able to 
rewrite them as a series of first order equations, as we have done above for the 
boson eigenvalue problem, by redefining the field derivative as a new variable, and 
with this definition serving as a a new equation. In each of the other equations the 
second derivative of the field is replaced by the first derivative of this new variable. 



and the first derivative of the field is replaced by the new variable itself [^. In 
general we then have to satisfy N equations of the form 
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%^ = Mx, y,, ...,yr,),i = l, 2, ..., N. (236) 

ax 

For simplicity consider A^ = 1. This can be solved by 

Vn+l =yn + Ax f\Xn, Vn) + 0{/\x^). (237) 

which moves the solution from x = Xn to x = Xn + Ax. This is the Euler method. 
To increase accuracy and stability we expand the second term in a Taylor's series 
about f{xn + Ax/2,y„ + Ay/2) (where Ay = Axf{xn,yn)) giving on writing the 
O(Ax^) term as well 

Ax^ 
Vn+i = yn + Ax f{xn + Ax/2, y„ + Ay/2) - -^f{xn + Ax/2, yn + Ay/2) (238) 

Ax^ 
+ -l^f'i^n + Ax/2, yn + Ay/2)/2 + 0{Axf . 

where the +|Ax^/'(x„Ax/2, yn + Ay/2) term comes from the O(Ax^) term of (237). 
This sets up a cancellation in the first order terms and we have a second order 
accurate method 

y„+i = y„ + Ax /(x„Ax/2, y^ + Ay/2) + 0{Ax)\ (239) 

This then is the spirit of the second-order RK scheme. The code uses a fourth order 
RK scheme. This uses derivatives at the initial point, at the half step, and at a 
full step, with weights of 1/6, 2/3 and 1/6 respectively as can again be confirmed 
by Taylor expansions. Guesses for the y value at the 1/2 step are made using the 
derivative at the original point and the derivative at the previously guessed point 
and these two are then averaged. For the full step the derivative used is that at a y 
value calculated using the second of the just described half step derivatives. 

2.3.2 BD Case 



Equilibrium configurations of boson stars in BD theory have been discussed in [ ^ . 
They used a modified version of the code we used to get equilibrium GR configura- 
tions. They added the BD terms to the equations as well as adding the BD scalar 
field. The equilibrium configurations in this case were those for which the boson 
field has the e*"^^* dependence in time as described before, while the metrics and 
the BD field are static. The value of the central BD field was not a free parameter 
but depended on its specified value at cxo. This was chosen to be zero to some tol- 
erance. The equations were integrated for different values of the central lapse value 
for a given central boson field value and guessed central BD field value. The value 
of the Brans-Dicke field at the edge of the grid was checked against its expected 
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value (~ ipoo + C /r) to the predetermined tolerance. If it did not match the value 
of the central BD field was changed and the process repeated for different central 
lapse values until a match was found. Since these configurations were quite similar 
to the GR case the central lapse values were close to the GR values. For all self- 
couplings the BD particle number was very slightly smaller than the GR case while 
the mass was slightly lower for small values of self-coupling and slightly larger for 
higher values of self-coupling. This was true for both the ground state and excited 
state configurations. 

2.3.3 Boson Halos 

The equilibrium configurations are obtained [BD equations with m = a = 0), from 
the equations 
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Here we use the dimensionless variables r = ujr, t = ujt and a = v47rG0. These 
equilibrium configurations are characterized by static metrics although the fields 
themselves again have an e*"^* dependence. Regularity at the center implies g{r = 
0) = 1. The equilibrium configurations are characterized by saddle points in the 
density p, which might have some significance in stability issues. These systems are 
characterized by a two parameter family of solutions. Firstly, we can vary o"(0), the 
central density and for each o"(0) we can vary A^(0). This is effectively changing the 
frequency and for each value of the frequency of the system there is a configuration 
for which the dispersive effects cancel the gravitational effects. We have noticed an 
interesting feature of the mass function, 



M{r) 




(244) 



where g^ is the radial metric. If plotted at a given value of r the profile is very similar 
to boson stars as shown in figure 4a. We also see the increase in mass linearly with 
radius for large r in the relativistic case. Figure 4b shows mass versus radius for 
different central densities and a given value of iV(0). 
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Figure 2 Mass profiles of ground state boson stars for different values of the self-coupling 
A are shown. The increase in mass with A is clear although the profiles are very similar. 
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Figure 3 Masses of zero-node, one-node and two-node boson stars without self-coupling 
are plotted as a function of central density. The maximum mass of one-node stars is 
1.356 M PI /m while the maximum mass for two-node stars is (as expected) greater at 
2.095 M PI /m. The profiles are deceptively similar to their ground state counterparts. 
Excited state stars are inherently unstable irrespective of the branch they lie on, unlike 
ground state stars that can be termed stable or unstable depending on whether they lie 
on the branch to the left of the maximum mass or to the right respectively. 
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Figure 4 (a) Left: The mass profile of equilibrium configurations of massless scalar fields, 
are interestingly similar to the profiles of massive scalar fields and neutron stars. The 
mass for a given radius increases to a maximum with central density before it decreases 
with further increase in central density However by calculating the particle number, we 
see no division into stable and unstable configurations at the peak. The particle number 
depends very sensitively on the initial value of the lapse function N, so that the clear cusp 
structure seen in the mass-particle-number diagram for boson stars cannot be derived 
here, (b) Right: The profiles are largely independent of radius although the mass itself 
increases with r. 
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Chapter 3 

Numerical Evolution of Bosonic Objects in 
Spherical Symmetry 

In the previous chapter we developed (from the action for a BD system with a matter 
Lagrangian that produces a scalar field equation for bosons, and for a spherically 
symmetric metric) the basic equations that govern the structure of an equilibrium 
boson star (see |2.1.1D . The GR case was a special case of these with the parameter 
a set to zero and discarding the BD field equation. (The BD field decouples from 
the physics). We now write down the evolution equations for the BD case. These of 
course include the time derivatives of the metric and the BD field. For convenience, 
so that we have no second time derivatives, we introduce new variables 



— df(rip) = — 



U^ = i-^dt{r^) = -dt{r^), (245) 



U^ = j^dt{r^) = ht{r^), (246) 



where ^ = ^/A7tG^(J). 

We thus get the field equations (for 0, \l/) 

dt{r^)=f3U^, (247) 

dtU^ = idr(3) drir ip)+/3 drdrir cp) - {dr/3) (r if) ^ 



r 



-Ngr 2ae^" 



- V,,^ V''^"^ + 2 e^" r (^^"f) 



(248) 



dtir^)=/3U^, (249) 



dtU^ = (drp) dJr ^) + /5 drdJr ^) - (dr(3) (r^) - 

r 



-2Ngr 



oadVi^^^) 



e-^^-aV.v^V> 



(250) 



where there are really two components for the second equation above (since we have 
a complex boson field). On each time slice we solve for the lapse A^ by using the rr 
component of the Einstein equation. However, instead of using the 00 component to 
get the radial metric on each time slice it is more convenient to use this component 
to record the time development of the radial metric through the Or equation (the 
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momentum constraint equation). The momentum constraint and Grr equations are 
respectively 

dtg = N [n^ drip + e^'^ ^ (n^ dr^ + n^ dr^% (251) 

A^ , o ,. Nr 
2r 



drN = ^ (/ - 1) + ^ ( {drvf + nj ^ 



+e 



2a 



{dr^) {dr^^) + n^ n; - - 2 g^ e^'^ V{^¥) 



(252) 



Since Too is the energy density the Goo equation is the Hamiltonian constraint 
equation. It is given by 



gr r^ 



S + (5^^)2 + e2a 



n*n], 



+ ((9,^) ((9^^"^) + 2 c/2 e2» \/(^^t) 



(253) 



This is used to monitor the code by keeping track of the accuracy with which it is 
satisfied throughout the evolution. Orders of 10"^ or better are typical at the end 
of our runs. 

3.1 Evolution Equations and Code 

The equations that describe the evolution of the boson star in GR are given by 
those in (249)-(252) with the parameter a set to zero, and again discarding the 
wave equation for the BD field. In order to finite-difference the derivatives, one 
must rely on well studied methods used in the study of linear differential equations, 
although we are in fact dealing with very non-linear systems. Our equations are of 



the hyperbolic form (wave equation) and we use a leapfrog method for evolution [51 



Leapfrog Method 

Consider a simple outgoing wave equation of the form 
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(254) 



Now take a one-sided differencing scheme for the time- derivative {n rep- 
resenting the timestep and j the position of a space coordinate) 
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+ 0(At). 
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For spatial derivatives consider a second-order representation 

u] = u]^^ - Axu'lj^n + {Axfu" + 0{{Axf) (256) 

u] = u]^, + Axu'\j,n + (Ax)\" + 0((Ax)=^) 

so that on subtracting the two equations 

u'kn= ^^\^^'-' +OiiAxn (257) 

Thus this 'Forward-Time Centered-Space' scheme gives 

^— L = -y ^+^ J-i . 258 

At \^ 2Ax J ^ ' 

This scheme is unstable as shown by Von Neumann stability analysis. 
Consider the coefficient v to be slowly varying enough to be considered 
constant. Now consider eigenmodes of the difference equation of the form 

C - 1 = -^ i sin(A;Ax). (259) 

The time dependence of the eigenmodes grows like .^" and since |.^| = 
a/I -f (■^)^ sin^(fcAx) > 1 this instability grows and renders the scheme 
unstable. 

So instead one uses a staggered leapfrog method using centered-time 
and centered-space derivatives. 

Now we get a quadratic equation for ^ with a solution |^| = 1 for any 
^^ < 1. Thus there is no growth of instability or amplitude dissipation 



51| , |37| . It is called a leapfrog scheme because the time derivative terms 
are on odd time slices while the time levels in the space derivative terms 
are on even time slices. 

When one has second-order time derivatives (as we do) it is more appro- 
priate to redefine the ffist time derivative in terms of another variable 
(tt = j^ Qt ) ^^^^ ^^ then given on the ri + ^ steps. 
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The two components of the equihbrium scalar field have components \E'i with 
time dependence as cos(ci;£;t) and \E'2 with time dependence as sm{u Et). Thus quite 
generally at t = we set the field \E'2 = 0. Our equilibrium code calculates the 
initial field configuration \l'i and the metric components as functions of r. This is 
the initial data for the evolution code. The first spatial derivatives of the metrics 
and fields are then determined by using f'Ax = a {f{x + 3Aa;) — f{x — 3Aa;)) + 
b (/(x + 2Aa;) — f{x — 2Aa;)) + c {f{x + Ax) — f{x — Ax)). The coefficients a, b 
and c are determined (a = 1/60,6 = —9/60, c = 45/60) from the conditions that 
the coefficient of the first derivative of the Taylor expansion be Ax and the third 
and fifth derivative coefficients be zero. This makes /' sixth order accurate since 
the even derivatives automatically vanish. In order to get sixth order accuracy in 
the second derivatives the form above is applied for /"(Ax)^ this time with every 
f{x — s) and f{x + s) added rather than subtracted so odd derivatives vanish. The 
coefficients are then determined (a = 2/180, b = -27/180 , c = 490/180). At the 
last but one grid point though one only has terms to second order accuracy. 

The metrics and fields have vanishing first derivatives at the origin. From the 
equations it is clear their parities are conserved. A pair of fictitious points corre- 
sponding to negative r are set up. The code evolves rip and tt which are antisym- 
metric about the origin. The derivative of rip at the origin is ip and that is how we 
determine it. The radial metric at the origin is unity but the lapse is not prescribed 
there. Instead we choose to fix the lapse at the outer boundary and integrate in- 
wards on each slice using a fourth order RK method. Rather than integrate the 
Hamiltonian constraint equation on each slice to determine g we monitor its evolu- 
tion from time step to time step. We only determine the boundary value of g from 
the Hamiltonian constraint equation. 

The outer region of the system is the low field region and the mass of the system 
can be determined using the value of the radial metric there. 



M = lim - 

1 — ►oo 2 



g'^ir) 



(261) 



The outer boundary condition on the fields is an outgoing wave condition. If in 
the asymptotic region the field is of the form ~ g«(k.r-a;t)^ then this results in the 
condition (from the KG equation) 



g2 - JY2 



m^ (262) 

living k = ^77(1 rn^ 
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which can be written (in terms of the fields and dimensionless units of the code) as 
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Ar2 
ii = -i3n' r^, (263) 

where (3 = —. This ehminates first order reflections off the grid. In order to remove 
higher-order reflections, we add a "sponge" in the outer region of the grid to absorb 
reflections. To the 11 equation (250) we add a damping term at the edge of the grid. 
This "sponge" is effectively hke adding a potential term that is proportional to A; + ti; 
for an incoming wave and proportional to u — k for an outgoing wave. Thus it forms 
a potential barrier against incoming radiation. This sponge is smoothly phased in 
from some point r(j) out to the edge of the grid. The width of the sponge is of the 
order of the wavelength of the scalar radiation. 

3.2 Nature of the Perturbations 

We study the dynamical properties of the boson stars by perturbing the equilibrium 
fleld distribution. The accretion or annihilation of scalar particles is simulated by the 
addition of fleld in the outer regions of the star or by decreasing it in denser regions 
of the star respectively. Another type of perturbation that has been studied is 
changing the ipi and il)2 of the equilibrium conflguration. This perturbation changes 
the kinetic energy density distribution. In either case the changes in the metric 
functions g and N are determined by the constraint equations and the polar-slicing 
condition. That is, we reintegrate the Hamiltonian constraint and lapse equations of 
(252) and (253) on the initial time slice. The magnitude and the length scale of the 
perturbations can be chosen arbitrarily. The perturbations are always spherically 
symmetric. 

3.3 Results and Conclusions for ID Boson Stars with Self- 
Coupling 

• Ground State Results 

An equilibrium conflguration of boson stars when perturbed reacts to these per- 
turbations in a myriad of ways. As discussed before, in the works of Pang and 
Lee H], ground state boson stars without self-coupling on the S branch have very 



speciflc quasinormal modes of oscillation under inflnitesimal perturbations. In |^ 
the mass function is plotted against radius for a given particle number. The curve is 
seen to have a "well" with a minimum at the exact equilibrium mass corresponding 
to that particle number. For any other given mass at that particle number there are 
oscillating regimes from a minimum to a maximum radius. In a numerical scenario 
we are in a position to test more general and realistic perturbations. While we do 
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study the effects of infinitesimal perturbations (by applying perturbations for which 
the mass change from equilibrium is less than 0.1%.) larger perturbations are also 
analyzed. The effects of larger perturbations cannot be studied analytically using 
perturbative methods. From the former we get quasinormal modes of oscillations, 
while from the latter we get scalar radiation that results in mass losses in these 
stars followed by their eventual settling to stable configurations of lower mass, their 
collapse to black holes, or their total dispersal. Once perturbed the metric equa- 
tions are reintegrated to get consistent metrics for those new initial configurations. 
Thereafter the system is allowed to evolve under the evolution equations. 

• S Branch Stars 

In the case of a free field (A = case) a perturbed S'-branch star oscillates with a 
definite frequency, losing mass through bursts of scalar radiation at each expansion. 
It finally settles to a new S'-branch configuration of lower mass. Here the effect of 
a self-interaction term on this behavior is examined. In Fig. 5 a typical example of 
the perturbed field configuration, and radial metric, of a star with A = 10 and a 
central density of o"(0) = 0.1 is shown. This star has been perturbed by accretion 
of scalar particles in a region of lower density. Its evolution is detailed below. 

Fig. 6 shows the radial metric as a function of distance for the same configuration 
at various times. The labels A, B, C, D correspond to times (in units of the inverse of 
the underlying scalar field frequency) t = 192, 306, 391, and 505 respectively. The 
positions of the peaks are labeled Rq, Ra, Rb, Rc, and Rd, where Rq is the posi- 
tion of the initial unperturbed peak. Here Rq = 7.95, Ra = 8.55, Rb = 6.6, Re = 
8.2 and Rr, = 6.65, where the length scale is in terms of the inverse mass of the 
boson. The oscillations are shown clearly in Fig. 7 which is a plot of the maximum 
value of the radial metric as a function of time. The point where this function is a 
maximum corresponds to the core of the star contracting to its minimum size in a 
cycle. Similarly, the maximum radial metric starts to decrease as the star expands. 
At each expansion the star loses mass through scalar radiation. The oscillations 
damp out in time as the star starts settling down to the new configuration. A plot 
of mass vs. time is shown in Fig. 8. The emission of scalar radiation decreases as the 
oscillations damp out, as can be seen from the figure. The slope of the curve steadily 
decreases as the star starts settling down to its new lower-mass configuration. The 
mass is measured at the inner edge of the sponge using (261). Exact details of the 
curves have some dependence on the sponge parameters but the basic results are 
the same. 

A characteristic of the boson star that might be important for its potential 
observation and identification is its fundamental oscillation frequency which can 
be determined from Fig. 9. We found / = 1/(199 A^(oo)) = 4.7 x 10"^. The 
oscillation frequencies for a large number of S- branch stars have been compiled in this 
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way. Fig. 9 shows a plot of the oscillation frequency versus mass for many slightly 
perturbed configurations (masses within 0.1% of the unperturbed mass.) As the 
mass increases, the frequency increases and then drops down as the transition point 
{dM/dcr{0) = 0) is approached, signalling the onset of instability. This is seen for 
both the non self-interacting as well as the self-interacting case. These quasinormal 
modes of oscillation characterize S'-branch stars. The point of transition from the 
5* to the U branch corresponds to a zero of the oscillation frequency. 

For a given mass higher A S'-branch stars have a lower oscillation frequency than 
similar mass lower A stars, unless one is near the transition point of the lower A 
configuration. This is not too surprising, since for a given mass the radius of the star 
increases with increasing A. We have seen this trend even for A values as high as 
1600 (see discussion on high A stars in section ( p.3.2| )). However since the maximum 
mass of higher A configurations is greater than that of lower A configurations their 
maximum oscillation frequency could be greater than that for lower A stars. (This 
can be understood as a size effect: on the S-branch, higher mass stars are smaller 
and have higher frequencies.) This can be seen in Fig. 9 . The maximum oscillation 
frequencies of A = 5 and A = 10 configurations are higher than that of the A = 
case. As the stars get much larger, though, the highest frequency starts to decrease. 
For example the maximum frequency for A = 30 stars is less than that of A = 15 
stars which is itself lower than that of the A = 10 case. A perturbation calculation 
for the high A case is presented in section ( p.3.2| ) to show the dependence of the 
quasinormal mode frequencies on A for high A configurations. The frequency is 
then proportional to the inverse of the square root of A. Thus as the stars become 
extremely large they oscillate less and less rapidly and it is no longer feasible to evolve 
them numerically in our code. (The time step used in the numerical simulation 
cannot be increased as it is determined by the intrinsic oscillation time scale of the 
scalar field, which is many orders of magnitude shorter than the oscillation time 
scale of the whole star for these cases.) 

The quasinormal mode curves are also useful in determining the evolution of 
strongly perturbed S-branch stars and the final configurations they could settle 
into. A perturbed star loses mass and settles to a final configuration corresponding 
to some position on the solid line in the figure. In Fig. 10 we single out the A = 10 
curve and plot the evolution of the S'-branch star discussed above. The points PI, 
P2, and P3 show the route to a new configuration. These points correspond to 
times t = 0, 1200, and 4800 respectively. By extrapolating this line to where it 
meets the A = 10 curve one could expect a final mass of about 0.76 Mpi/m . 

• [/-branch Stars 

Ground state ^-branch stars (with or without self-coupling) can migrate to the S 
branch under certain perturbations. In general they are unstable and cannot settle 
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into new configurations on the same branch. In [^], with the aid of catastrophe 
theory, the instability of these unstable stars is also discussed. Their conclusion 
regarding unstable boson stars is that they disperse or collapse. However we also 
demonstrate the presence of migrations of these stars to configurations on the S 
branch. This is due to our perturbations being more general, allowing for changes 
in mass and particle number, whereas they keep them strictly constant. 

Figs. 11-13 show a migrating A = 30 star whose unperturbed overall field density 
a has been decreased by about 10%. Fig. 11 shows the behavior of the radial metric 
in time as a function of radius. It oscillates about the final ^-branch configuration 
that it will settle into. This final state is shown on the plot as a dark line. Fig. 12 
shows the maximum radial metric as a function of time. The star initially expands 
rapidly as it moves to the S'-branch. This can be seen from the sharp drop in 
the radial metric. Once it moves to the stable branch, it oscillates about the new 
configuration that it is going to settle into. Fig. 13 shows the mass of the star as a 
function of time. It loses mass at each expansion losing less and less mass at each 
subsequent expansion and the curve gets smoother and smoother as it settles to its 
final state. 

Fig. 14, which shows the oscillation frequency as a function of mass for A = 30, 
can be used to predict the end point of migration. Points Ql, Q2, Q3, and QA show 
the migration of this star. These correspond to times of 500, 1000, 2000, and 3500. 
The oscillation is clearly damping out. The final configuration it is expected to 
settle down to is shown as a dot and corresponds to a stable star of central density 
cr(0) = 0.0817 with a mass of 1.037 Mpi/m. This example is typical of a number of 
simulations of [/-branch ground state configurations with self coupling. 

For higher central density stars on the ^/-branch, the mass versus central density 
curve has a second, gentler peak, similar to the white dwarf neutron star situation. 
One might suspect that this corresponds to another stable and unstable branch 
respectively. However we find that configurations on both sides of the peak are 
unstable. These configurations always disperse upon perturbation, consistent with 
the fact that they have M > N m, where M is the mass of the star and N m is the 
number of bosons multiplied by the mass of an individual boson. 
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Comparison of Unperturbed and Perturbed Initial Profiles 
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Figure 5 Comparison of a strongly perturbed ground state A = 10, S-branch star [M = 
0.781 Mli/m, o-(O) = 0.1] to the unperturbed configuration [M = 0.722 MlJm]. The sohd 
hnes correspond to the unperturbed configuration and the dashed ones to the perturbed 
star. The perturbation shown corresponds to the addition of scalar field a at t = 0. 
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Figure 6 The evolution of the radial metric gr 



g^ for the configuration shown in 



the previous figure. The initial perturbed configuration labeled t = and unperturbed 
configuration are shown. A, B, C, and D correspond to times t = 192, 306, 391, and 
505 respectively. The positions of the peaks for these times are labeled Ra, Rb, Rc, and 
Rd- 
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Figure 7 The peak value of the perturbed radial metric is plotted over a long time 
(6000 h/mc). The points labeled A, B, C, and D correspond to the same labels in the 
previous figure. The oscillations decay in time. 
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Figure 8 The total mass of the star is plotted as a function of time. The mass loss 
through scalar radiation decreases in time as the oscillations start damping out. 
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Figure 9 The oscillation frequencies of different ground state boson star configurations 
are plotted as functions of mass, for A = 0, 5, 10, 15, 30, 100, and 200. The curves are 
obtained by slightly perturbing (perturbed mass within 0.1% of the unperturbed mass) 
S-branch stars. They reach a peak and then drop down at the approach of the maximum 
mass allowed for a given A, signalling a transition from stability to instability. 
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Quasinormal Modes and Stable Star Migration 
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Figure 10 The highly perturbed A = 10 S-branch star has an osciUation frequency below 
the A = 10 solid line. Its movement towards the solid line is shown through points PI, 
P2, and P3 corresponding to times 0, 1200, and 4800 (fi/mc) respectively. 
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Figure 11 The radial metric g^ = §„ of a perturbed U -branch ground state star is 
shown at various times. The curves 1,2, 3, a, b, c, and d correspond to times t = 
10, 20, 30, 340, 440, 540, and 640. The unperturbed star has a central density a{0) = 
0.23 and a self coupling parameter A = 30. The initial equilibrium metric configuration 
has also been shown. The initial field density of this star has been lowered by about 10%. 
The t = curve corresponds to the initial perturbed radial metric. In the asymptotic 
region g^ is not oscillating but monotonically decreasing due to the mass loss. 
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Figure 12 The maximum value of the radial metric is plotted as a function of time. The 
initial sharp drop in the radial metric is due to the expansion of the star as it proceeds to 
the stable branch. There it oscillates about the new stable configuration that it will settle 
into. This corresponds to a star of mass M = 1.037 Mp^/m (whose metric configuration 
has been shown in the previous figure as a dark line). The points a and c correspond 
to two minima in the peak of grr which occur when the core of the star reaches its local 
maximum size. Likewise the maxima in the peak of grr at b and d correspond to the core 
of the star reaching its local minimum size. 
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Figure 13 The mass is plotted against time. The mass loss through scalar radiation 
at each expansion of the core (corresponding to the maximum radial metric reaching a 
minimum) decreases in time as the oscillations damp out. 
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Quasinormal Modes and Unstable Star Migration 
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Figure 14 The migration of the U -branch star considered in the previous figure is shown 
after the star has moved to the S-branch. Points Ql, Q2, Q3, and QA correspond to times 
500, 1000, 2000, and 3500 (h/mc) respectively. Here too the mass loss decreases in time 
and the star finally settles to a stable configuration. 
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3.3.1 Excited States 

The similarity of mass profiles of excited boson stars to their ground state coun- 
terparts, (Fig. 3) might lead one to expect stable and unstable configurations to 
the left and right of the maximum mass respectively in analogy with ground state 
configurations. However, our numerical studies show that the excited boson star 
configurations on both sides of the peak are inherently unstable except that the time 
scales for instability are different. If they cannot lose enough mass to go to the 
ground state they either become black holes or totally disperse. This occurs even 
if no explicit perturbations are put in the numerical evolution other than those in- 
troduced by the finite differencing error in the numerical integration. We have also 
carried out perturbations that correspond to scalar particles falling onto the star, 
or those that decrease the scalar field strength at the center point (corresponding 
to scalar particles decaying through some channel). The instability shows up in all 
cases studied. We note that this instability is NOT in contradiction with the result 
of [Q, which concluded that S branch excited states are stable under infinitesimal 
perturbations that strictly conserve M and A^, (where M is the mass of the star, 
A^ the number of bosons and m the boson mass.) Our result of instability under 
generic perturbation is consistent with the studies of A = stars under infinitesimal 
perturbation [Q. The presence of a self-coupling term increases the time scale of 
instability but the essential pattern remains the same. Fig. 15 shows a perturbed 
1-node star (whose mass has been reduced by about 8% to 0.9 Mpi/m by a pertur- 
bation) making a transition to the ground state although the mass is greater than 
the maximum ground state mass of 0.633Mpi/m. A substantial amount of scalar 
radiation is emitted in the dynamical evolution, which brings the mass below the 
critical value. The evolution of the radial metric function is shown. That there are 
two peaks at t = is indicative of a first excited state. One of the peaks disap- 
pears gradually as the star goes to the ground state branch. The star then oscillates 
about the ground state configuration that it will finally settle into. In Fig. 16 a 
three-node configuration with a total mass of of 0.92Mpi/m is shown going to the 
ground state branch, after radiation by scalar waves carries off the excess mass and 
kinetic energy. The plot shows the density function against the radius and time of 
evolution. By the density function we mean density p multiplied by r^, that is, the 
mass per dr at radius r. The density function, pr^, has n + 1 maxima for an n node 
star and hence here we have four sets of lines initially. At the end of the simulation, 
we see that it settles down to a ground state configuration with small oscillations 
with ever decreasing amplitude. For these simulations of low central density stars 
we put in an explicit perturbation to the equilibrium configuration since otherwise 
the instability time scales are extremely long. 

For stars with higher central density, there is a critical density above which 
the stars cannot lose enough mass to go to the ground state but collapse to black 
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holes. In our numerical simulation for one node stars this critical density is o"(0) = 
(72 = 0.048. As the central density is increased the kinetic energy of the highly 
compressed initial equilibrium configuration is increased and the star first expands 
before the eventual collapse to a black hole. As the central density further increases 
towards the M = Nm point the expansion phase becomes longer. In Fig. 17 we 
show the density function (pr^) against radius at various times for 4 configurations 
(an S'-branch configuration with o"(0) = 0.1, and three [/-branch configurations of 
central densities o"(0) = 0.3, o"(0) = 0.4, and cr(0) = 0.5). The initial configurations 
are the equilibrium ones without any explicit perturbation (except those introduced 
by the discretization used in the numerical simulations). The first frame shows 
the ^'-branch star radiating a little as it almost makes a transition to the ground 
state. However it cannot sustain this state for long and it collapses to a black hole. 
The decay time decreases in the case of a [/-branch stars of central density 0.3 
and 0.4. However for the o"(0) = 0.5 star, the star is clearly getting more diffuse 
than the previous ones. It goes through an expansion phase initially although it 
finally collapses to a black hole. Configurations with o"(0) > 0.54 have M > N m. 
They do not collapse to black holes but disperse to infinity. Next we turn to a 
study of instability time scales. In simulations where a black hole will form, the 
imminent development of an apparent horizon leads to a rapid collapse of the lapse 
due to the polar slicing used in the evolution. We take the time of collapse of the 
lapse at the origin to ~ 10^^ of its initial central value to be the approximate time 
for formation of the black hole. Fig. 18 shows this time scale for a 1-node star 
without self coupling. We plot the decay time scales of first excited state stars as a 
function of central field density. Again no explicit perturbation is applied in these 
evolutions other than the discretization error in the simulation. In order to make a 
fair comparison of the time scale due to such a perturbation we cover the radius of 
the star in all cases by the same number of grid points. The maximum ground state 
mass for stars without self coupling is around 0.633 M|,;/m. This corresponds to a 
central density of o"i = 0.021 for a one-node star. We mentioned earlier that stars 
with central densities below o"2 = 0.048 (mass of 0.91Mpi/m) lose enough mass and 
move to the ground state. Beyond that and up to the point M = N m they collapse 
to black holes, while stars with M > Nm, corresponding to a central field density 
of cr(0) > 0.541, disperse to infinity. The time by which collapse takes place to a 
black hole decreases with increasing central density along the S branch. This trend 
continues for a while into the U branch (starting at o" = 0.25 until a ~ 0.4), but as 
one approaches the M = Nm point the stars lose a significant amount of matter to 
infinity before they collapse to black holes and evolve on a longer time scale. For 
example, a star of central density 0.5 has an initial radius of r ~ 9. (Remember 
that the radius of the star is defined as the radius which contains 95% of the mass of 
the star.) Its radius increases to as much as 115, more than an order of magnitude. 
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before it starts collapsing. The dispersion time scales of a couple of typical examples 
of stars for M > Nm which disperse to infinity (instead of collapsing to a black hole) 
are also shown on the figure. To give a sense of the instability time scales of these 
configurations we take the time scale to be the time by which these stars disperse to 
ten times their original radii. This time drops drastically for stars with M ^ Nm. 

Next we turn to the case of A 7^ 0. Fig. 19 demonstrates black hole formation 
for a A = 30 star in the first excited state, with a central density cr(0) = 0.1. This 
star had an initial radius of about 20.7 where the radius is again defined as that 
containing 95% of the mass. The lapse finally collapses towards zero indicating that 
an apparent horizon is about to form. The time scale of collapse to the black hole 
was around 1985 compared to a time scale of less than 800 for a o"(0) = 0.122 star of 
similar radius without self-coupling which is shown in Fig. 20. This is to be expected 
as the A term represents a repulsive force. The collapse time is again defined by the 
lapse collapsing to 10~^ of its original value at r = 0. 

We now turn to the evolution of a highly excited state. In Fig. 21 we show the 
initial field configuration of a star containing five nodes. For a five-node star the 
density has a central maximum, and then five local maxima, each subsequent one 
smaller than the one preceding it. This star is then evolved without any perturbation 
except those introduced by the discretization error of the numerical evolution. In 
Fig. 22 we show a contour plot of the evolution of the density function in time: pr^ 
has n + 1 maxima for an n node star and hence here we initially have six sets of 
lines centered at dimensionless r = 5, 13, 35, 52, and 75 respectively. This star has 
central density cr(0) = 0.075 and it collapses to a black hole after a long evolution. 
In the process we see intermediate states with fewer numbers of nodes. 

For comparison, in Fig. 23 we show the contour plot for the five-node star up to 
a time of t = 1000, at which time it has decayed into a four-node state, against the 
equilibrium density function of a four-node star with central density a = 0.06. Very 
clearly the maxima are at similar radii (although the size of the peaks are somewhat 
different). 

This feature, that of nodes disappearing and the star cascading through lower 
excited states, is characteristic of the decay of higher excited states of boson stars. 
Although in this particular case and at this particular time the decaying star is close 
to a specific lower excited state, in general the decaying star is (roughly speaking) 
a combination of lower excited configurations. This is similar to the decay of atoms 
in excited states. However, we note that for the "gravitational atom" [^ there is 
no exact superposition, due to the intrinsic non-linearity of the system. 

Whenever we report black hole formation, we are of course talking about immi- 
nent formation since our highly singularity-avoiding polar slicing condition causes 
the lapse to collapse as soon as an apparent horizon forms. However taking the data 
just before such formation, and evolving it with a 3D code with maximal-slicing has 
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allowed us to confirm that a black hole indeed forms. 

3.3.2 Very High A Configurations 

While calculating the eigenvalue frequency for the equilibrium boson star it is found 
that it gets increasingly difficult to calculate the eigenvalue as the value of A gets 
large, because the equations become very stiff. There are two scales to the problem: 
a scale of slow variation of the field inside a certain radius (related to A) followed by 
rapid decay outside it. This produces an effective surface layer to the star making 
it similar to neutron stars. It turns out that the large A limit can be treated using 
a set of approximate equations that are exact in the A = oo limit [6]. By making 
the following change of variables: 
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where the primes refer to differentiation with respect to f. In the limit of A 

one can keep terms to leading order in -^ in (265) -(267). In particular (265) 

reduces to an algebraic constraint 
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However, we note that this is valid only for a ground state configuration. For a 
state with nodes it would not be reasonable to neglect derivative terms compared 
to terms proportional to o", because it is zero at every node. 

To get an estimate of the accuracy of the high A approximation we compare 
the solution using the approximate equation for high A to the brute force numerical 
solution of the complete set of equations, for a A = 800, a = 0.05 boson star in 
the ground state (see Fig. 24). The agreement of the fields is quite good until 
the outer region where the approximate equations cause the field to abruptly fall 
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to zero. Comparisons of the radial metric and the lapse are also shown. In the 
A = oo limit the approximate equations are exact and the star really has an outer 
surface reminiscent of a neutron star. In fact, an effective equation of state can be 
written [6]. In Fig. 25 the mass and particle number versus central density a for high 
A stars is shown. One expects (as in the case of other ground state configurations) 
that the configurations with M > Np-m disperse when perturbed, while those on 
the U branch with M < Npirn would be unstable and, if unable to migrate to the 
S branch under perturbations, would form black holes. (Here we use the symbol 
Np for the particle number and not A^ so as not to confuse it with the lapse as 
both these functions figure prominently in the analysis that follows.) The particle 
number is calculated from the current J° and is given by 



Np = An r^j-iua'dr. (269) 

Here a has been replaced hj uj a and d^r by Anr'^dr for the spherically symmetric 
case. In terms of the bar coordinates we see that Np = yKNp. 

Next we turn to the determination of the quasi-normal frequency (QNM) of 
the high A stars. In principle one could determine the QNM using the dynamical 
studies as performed for the A = case. However the procedure is extremely 
computationally expensive. The scalar field has an inherent oscillation period of 
about 27r and the evolution time steps must be small enough to resolve it. However 
the code must run long enough to see a few metric oscillations in order to determine 
the quasinormal mode. As A gets large, the sizes of the stars also get large leading 
to a lower frequency of oscillation. In order to determine the QNM we use instead 
the following perturbation analysis based on [8,11] but using our own notation for 
lapse, fields, time, radius and self-coupling as defined in section II. 

We write the perturbed fields as: 

a={a^+ia2)e'^\ g = go + 6g, N = No + 6N (270) 

where 

ai = ao{r){l + 5ai{r,t)), a2 = ao{r) 6a2ir,t). (271) 

the Klein-Gordon equation can be written as 
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Calculating 6&2 from (279) and (282), and substituting into (274) gives, along 
with equation (283), 
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Adding 6Tq to 5T2^, and substituting in equation (283) and its derivative, we get 
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Using the expression for the particle number A''^ = /q°° d^xJ^^ where J^ 
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In terms of bar coordinates, defined in the beginning of this section, equation 
(284) becomes 
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Similarly equation (286) becomes 

90 -2^ I 1 



SNp = 4n r dr r^ ^4^1 -^-^ {go Sg' - go' 6g) (289) 

Jo No I go^ r a^ 




^ '^5-a' 



Ai-5 ao 




No' go' A \ao 

Now using SNp = (charge conservation), which is appropriate for large A where 
6Np is given by (289), we get on putting this in (287): 

Ai-5 ^ Ai-5 ^ yf No go J AO-5 A^g^ v / 

,^(_,„^(l,|),l^(i,^_w^^ 



AO-5^" M2A "M 2// Ai-5 do 



This equation suggests that 5g goes like 1/a/A- Writing 5g as —'x'Sg and (5o" as 
— X^<5o" we then see from (288) and (290) that the quasinormal mode frequency x for 
high A configurations must go like 1/a/A for a given a. We have numerically evolved 
stars with the same a and compared the QNM frequencies obtained to the inverse 
ratios of the square root of their A values for A = 800, A = 1200, and A = 1600 
confirming the above analysis. In Table I we show the comparison between the 
perturbation analysis and the numerical result. We see that the perturbation result 
gets more accurate for increasing A. We note that configurations that have the size 
of a neutron star would have to have A of order 10^^ with QNM frequency of order 

X = 10-^^ 
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Figure 15 The transition of a first excited state star to the ground state is shown. Here 
the radial metric is plotted against radius for various times starting from t = and 
then for ta = 250, t;, = 500, t^ = 1245, td = 4000, t^ = 5000, tf = 5505, tg = 6000, 
and th = 6370. The initial unperturbed and perturbed configurations are shown. (The 
perturbed configuration has the lower second peak.) The initial mass of the star was 
M = 0.901 Mpi/m after perturbation. The radial metric initially has two peaks indicative 
of a one-node configuration. As the star evolves and goes to the ground state one peak 
disappears. This can be seen in the curves from tc-t^. Once on the ground-state branch 
it oscillates and finally settles into a stable ground-state configuration. 
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Figure 16 The transition of a three-node configuration of mass 0.91 Mpi/m. This star 
loses enough mass during the course of its evolution to move to the ground state. 



94 









""»■ 


! 1—1=1025 














1-150 
























1-650 
























■ -1=1 025 




fe 


,; 


,, 


^£1==-=; 



Ij 




am-llA 


'1 










1 1 






1-100 




I 






1-150 






A 




1-1 BO 

— t-190 

1-210 






r'\ 


\ 






k 


' ^M 


/« 


\^'^>.-^-^ 




- 1-200 

- 1-300 

- 1-500 

- 1-700 
-■ 1-900 



Figure 17 A comparison of the manner of black hole formation of four excited state 
conhgurations. The first frame is an S branch star of central density (j(0) = 0.1, that 
tries to go to the ground state but fails to. As the central density increases, decays to 
black holes occur on a shorter time scale. A plot of the collapse of a U -branch star of 
central density o"(0) = 0.3 is shown in the next frame. Stars get more diffuse as one moves 
farther along the U branch. A star of central density 0.4 (shown in the third frame) has a 
decay time close to the previous one. Decay times then start to increase. A star of central 
density 0.5 (close to the point M = N m with cj(0) ~ 0.541J disperses to over ten times 
its radius before collapsing to a black hole. 
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Figure 18 The decay time to black holes is plotted as a function of central density, 
for one-node conhgurations (first excited states) of boson stars without self coupling. 
Perturbations are due only to the finite differencing effects of the numerical scheme. To 
make the comparisons meaningful the 95% mass radius (which is our definition for radius 
of a star) of every configuration considered was covered by the same number of grid points. 
Configurations for which the central density o"(0) < (T2 move to the ground state. The 
central density di corresponds to a mass M = 0.633 Mp^/m, which is the maximum mass 
of a ground state boson star. The decay time decreases with increasing cj(0) and this 
continues even for a(0) > 0.25 which is the transition point from the S-branch to the 
U -branch. The decay time then starts to increase as one approaches the M = N m point 
corresponding to a central density a{0) = 0.541, beyond which the stars disperse to infinity 
rather than becoming black holes. The dispersion times of two such stars to ten times 
their original radius (95% mass radius) are also shown on the figure. 



1.0 






0.5 



0.0 




■t^=1985 



-tb=1950 







20 



40 



60 



Figure 19 The evolution of the metric function N'^ = —gtt for a A = 30, cj(0) = 0.1 
boson star in the first excited state, without any exphcit added perturbation, is shown. 
The configuration lies on the S-branch and has an initial mass M = \.743>MpJm. The 
various time slices correspond to times t = 0, ta = 1060, t^ = 1950, and tc = 1985. The 
lapse function collapses as an apparent horizon is approached, signaling the formation of 
the black hole. (This is indicative of an inherent instability of the excited states). 
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Figure 20 The evolution of the metric function N'^ = —gtt for a A = 0, o"(0) = 0.122 
boson star in the first excited state, without any exphcit perturbation added, is shown. 
The configuration hes on the S-branch and has an initial mass M = 1.23 Mpi/m. The 
various time slices correspond to times t = 0, ta = 450, ti, = 720, tc = 750, and t^ = 770. 
This star has a radius of 20.7 which is about the same as that of the configuration shown 
in the previous figure. Again, the lapse function collapses when an apparent horizon is 
approached, as a black hole is being formed. The time scale of collapse is much less than 
for the A = 30 case in the previous figure. 
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Figure 21 The initial Eeld configuration of a five-node star. Tlie field has five nodes or 
extrema. The absolute value of each extremum is clearly smaller than the one preceding 
it. 
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Figure 22 A contour plot of a perturbed five-node star that decays into a black hole 
showing p r^ as a function of distance (vertical axis) and of time (horizontal axis) is shown. 
The density p is highest at the origin and has five other local maxima, each smaller than 
the previous one. The value of the maxima at the end are very small compared to the 
earlier ones, and to enhance the features p r^ rather than just p is shown in the plots. 
Each set of lines represents the maxima of p r^ and the number of lines in a set gives 
an indication of the height of the maximum. This particular conhguration has a mass 
M = 3.07 Mpi/m. This star cascades through an intermediate four-node state (around 
t = lOOOj before proceeding to form a black hole. Cascades are characteristic of excited 
boson star decays, similar to atoms in excited states going through intermediate states 
when decaying to the ground state. 



5 Node Star Evolution and Equilibrium 4 Node 
Configuration Cascaded To 




Figure 23 The equilibrium density function of a four-node star of central density a = 0.06 
(right frame) is placed alongside the contour plot of the Eve-node star described in the 
previous figure up to a time of t = 1000 when it has decayed into a four-node state (left 
frame). This plot shows that the transition of the five-node star is to a perturbed four-node 
state close to the one shown before it continues its evolution to a black hole. 
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Table 1 The ratio of the QNM frequency for A = 1600 to the QNM frequency for a given 
A is compared to A'^'^/40 (which is the predicted ratio for large A) for A = 1200, 800, 
and 1600. As expected, the higher A values match better. The initial central density is 
a(0) = 0.4. 
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Comparison of Exact and Approximate Functions 
For a High Lambda Configuration 
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Figure 24 The equilibrium profiles of a A = 800 star with central density a = 0.05 as 
derived from the high-A approximate equations and the exact equations are compared. 
A Schwarzschild exterior is attached to the approximate solution after the field vanishes. 
The three plots show the field a, §„, and goo respectively. Clearly, the approximation 
matches the exact solution very well. 
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Figure 25 The mass of a high A star generated from the approximate equations is plotted 
as a function of a (a/A^). It shows the same basic structure as the profiles generated 
for low A using the exact equations. The peak is at about 0.22A^''^ Mpjm,, which means 
that to achieve O.IMq we would have to take A of order 10^^. This would be a very large 
star to evolve numerically. Also plotted is N m (N is the particle number and m the mass 
of a boson). The crossing point of the two curves represents transition from negative to 
positive binding energy. 
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3.4 Boson Stars: Formation and Stability in BD Theory 



We now describe our work on the study of boson stars in the BD theory |^ . We 
hope to extend this work to general Scalar-Tensor theories in the future. 

The equations and set up have already been described earlier. We use the evo- 
lution code described in the GR section of this chapter to study the dynamical 



evolution of this system. The equilibrium configurations are described in [52|. A 
modified version of the equilibrium code used in GR was used to find the equilibrium 
configurations. We then studied the time evolution of these configurations using our 
evolution code. 

3.4.1 Boundary Conditions 

For the boundary conditions, regularity dictates that the radial metric is equal to 
1 at the origin. The boson field and the BD field are both specified at the origin. 
The boson field goes to zero at cxo and the BD field goes to a constant which is 
fixed during the evolution. This constant does not enter into any of the evolution 
equations as all the terms in the set of equations are derivatives in the BD field. 

The inner boundary at the origin requires that the derivatives of the metric and 
all the fields vanish at this point. This is implemented by extending the range of 
r to include some negative values. The metric components g, N, the boson field 
and the BD field are required to be symmetric about r = 0. While the boson field 
falls off exponentially with r in the asymptotic region, the BD field is more slowly 
varying with a 1/r dependence. 

For the BD field, ip, the wave equation can be written as 

c^(r,t) = c^oo + -F(t-0 (291) 

which is an exact condition for 1-D waves. Differentiating, one gets 

= 0. (292) 



1 dip d(p (if- v^oo) 



c dt dr r 



outer edge 



Note that this technique used by Novak ^Tj in a study of stellar collapse 



For the boson field, \1/, an outgoing boundary condition to order 1/r is 

/32A;2 = a;|-iV2m2e2". (293) 

Therefore at the outermost grid point we require that 

dtdt-^ = -13 dtdr^ - — ^ ^\ (294) 
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where \1/ = r \l/ is the field variable whose evolution equations are set up in the code. 
We also use a sponge to remove second order refiections. It reduces the momen- 
tum of the boson field artificially, and is irrelevant for the massless BD field. It is 
a damping term that effectively simulates a potential term in the wave equation. 
This potential is large for incoming waves (proportional to k + uje) and small for 
outgoing waves (proportional to fc — uje)- Therefore, we add an additional term in 
the evolution equation for 11^ (250) 

^^^(n^ + 9.^) (295) 

for Tend — D < r < Vend wlicrc Tend IS the r value of the outermost grid point and 
D is an adjustable parameter representing the width of the sponge. D is typically 
chosen to be a few times the wavelength of the scalar radiation moving out. 



The equilibrium sequences have similar profiles to GR [§2| with of course the 
additional BD profile. The perturbations are described in the GR subsection. While 
the mass loss in GR must solely be due to scalar field radiation there is an additional 
possibility here of gravitational field radiation through oscillations of the BD field 
even though we have spherical symmetry. 

3.4.2 Ground State Results 

• S-branch evolutions 

After confirming that our code was stable and the equilibrium profile of the star 
was undisturbed for several thousand time steps, we let the system evolve under 
infinitesimal perturbations (less than 0.1% mass change from equilibrium). The 
system started oscillating at its fundamental quasinormal mode frequency. The BD 
field also acquired this oscillation frequency. 

In GR, a QNM frequency increases as $c become larger (radius become smaller) 
up to a point before starting to decrease rapidly to zero as it approaches the density 
corresponding to maximum mass signalling the onset of instability |3^, Q. We 
found the same feature in BD. 

In addition, the system takes on the proper underlying frequency, which origi- 
nates from the time dependence of the boson field "^ ~ e**. This frequency corre- 
sponds to period 27r in t in our units, and leads to metric and BD field oscillations 
with period vr, from the structure of the equations. To show this feature we have 
enhanced the BD field oscillation (Fig. 26) in Fig. 27. The time interval between 
the ticks on the horizontal axis is vr. 

Under large perturbations, a stable boson star in GR expands and contracts, 
losing mass at each expansion. The oscillations damp out in time and the system 
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finally settles down into a new configuration on the S'-brancli. These features are 
now also observed in BD theory. 

We show the effects of a large perturbation on a stable BD boson star described 
above. The example we present is the case of initial data with $c = 0.2, mass M = 
0.54:0G^:/m after perturbation (about 13% lower in mass compared to unperturbed 
equilibrium configuration of M = 0.622G^,/m). 

The maximum radial metric grr and the central BD field as a function of time are 
shown in Fig. 28 and Fig. 29, respectively. In both figures, we have plotted the case of 
^BD = 600, and ujbd = 60 as well as GR. The increase of the maximum g^r indicates 
the star is contracting, reaching its maximum value at the end of the contraction 
in a cycle. Then, as the star expands, the maximum g^r decreases, reaching its 
minimum at the end of the expansion. These processes repeat themselves with the 
oscillations damping out in time as the star settles to a new stable configuration 
with maximum radial metric of smaller value than it started with (lower mass). The 
lower value of the BD parameter shows a phase shift in comparison to GR, which 
might be suggestive of a different rate of approach to the final configuration. We see 
the same dynamical behavior in the BD field (Fig. 29). It has the same oscillation 
frequency as that of the metric. The oscillations damp out in time and the BD field 
settles to a value closer to zero than it started at (lower final mass). 

The system loses mass through radiation during its evolution. A comparison of 
the mass as a function of time for the BD case {ujbd = 600 and 60) as well as GR 
shows little difference, indicating that the radiation is mostly scalar field radiation 
from the boson field and not scalar gravitational radiation due to the BD field. This 
is despite the BD field oscillating in the BD case and being zero in the GR case. 

The amount of mass radiated progressively decreases, as can be seen in the 
luminosity profile {—dM/dt versus time) shown in Fig. 30. Here again a comparison 
between GR and the two BD parameters is shown. Again we see the phase shift 
for the lower BD parameter (further from GR). The system finally settles to a new 
state of smaller mass. 

• Evolution of [/-branch stars and excited stars 

Boson stars on the U branch and excited states are inherently unstable in GR 
P3| , 0. Under perturbations that reduce the mass, boson stars on the U branch 
can migrate to the stable branch. In this section, we will see that these features also 
exist in BD theory. Configurations can be perturbed in such a way so as to decrease 
their mass enough that they migrate to new configurations on the S'-branch. On 
the other hand, in GR, if they do not lose mass and migrate, then the boson stars 
of configurations with M < Npirn collapse to black holes. Stars with M > Npin 
disperse and radiate out to infinity. We also see these features in BD theory. 
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Our first dynamical example from [/-branch boson stars is a migration process. 
As in GR, we have also seen migrations of these stars to the stable branch when we 
remove enough scalar field smoothly from some region of the star so as to decrease 
the mass by about 10%. In particular, we show the migration of a star of central 
boson field $c = 0.35 with unperturbed mass 0.625 G^/m. After perturbation, its 
mass is reduced to 0.558 G^/m. In Fig. 31, we show the maximum value of metric 
Qrr and central value of BD field (^{r = 0) versus time. The initial sharp drops 
in both lines occur as the star rapidly expands and moves to the stable branch. 
After that it oscillates and finds a new configuration to settle into. The damping of 
oscillations as it settles down is clearly seen in the figure. We see also the BD field 
oscillations damping out as the star gets closer and closer to its final state. 

The ratio of mass at time t to the initial mass for BD with parameter ubd = 60, 
600 and the GR case is shown in Fig. 32. The flattening of the curve at later time is 
indicative of the star settling down to a new configuration. Although convergence 
towards GR with increasing ubd is clearly indicated, there is no significant difference 
between the three cases. The amount of total mass extraction from the system is 
slightly suppressed if we evolve in the BD theory. By the time of 7500 shown in the 
plot, we see that the mass of the star is about 0.045 G^/m, which corresponds to 
an equilibrium configuration with $c = 0.06, while our central density $c is about 
0.061, meaning that the star is quite close to its final configuration. 

In contrast to the previous example, if we add a small mass to [/-branch stars, 
we can see the formation of a black hole in its evolution. In Fig. 33, we plotted 
an example of such an evolution, suggesting formation of a black hole. The initial 
data of this plot is boson star of central boson field $£ = 0.35 (the same value with 
the previous migration case) and the system is perturbed very slightly so that the 
perturbed mass is less than 0.5% greater than its initial mass. The sudden collapse 
of the lapse function is indicative of the imminent formation of an apparent horizon. 



(Due to the polar-slicing condition in our code |38| we will never see the horizon 
actually form). In addition to this the radial metric starts to grow rapidly and the 
code is no longer capable of handling the resulting sharp gradients. As an indicator 
of the suddenness of the process, we see that in the configuration shown the lapse 
has fallen to a value of about 0.003 by a time of 60 after being at 0.230 at a time of 
55 and 0.5 at a time of 50 (the latter two points are not shown in the plot). 

There is almost no loss in mass in this system and the time of collapse is quite 



close to that in GR. In the GR case, we confirmed the formation of black hole p8 



very shortly after this point (f» 3 M where M is the mass of the system) by switching 
this data into 3-dimensional code and evolving the system. Therefore we expect 
almost the same behavior in the BD case. 

The black hole formation in the BD theory has been investigated by Scheel, 
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Shapiro and Teukolsky for dust collapse 15^. They found that the dynamical be- 



havior of the apparent horizon is quite different in the physical Brans-Dicke frame, 
but the same as GR in the Einstein frame. Their results therefore also support our 
discussion of the formation of black hole in the BD theory. 

3.4.3 Excited States 

Excited states of boson stars in general are not stable in GR. They form black holes 
if they cannot lose enough mass to go to the ground state. We confirm the same 
features for the BD case. In Fig. 34, we plot the metric Qrr of dynamical transition 
from an excited state with 3 nodes to a ground state boson star configuration. The 
initial configuration has four metric maxima and the final has one. After it goes to 
the ground-state branch, it oscillates and compactifies to form a new configuration. 
We show the oscillations of the star from a time of 27300 to 28400 in Fig. 35. The 
95% mass radius at this stage is about 100 and the star has still to reach its final 
state. 

3.4.4 Formation of Boson Stars 

In the previous sections, we have analyzed boson stars in BD theory, starting with 
equilibrium or perturbed equilibrium configurations. However, we have not dis- 
cussed whether or how such an equilibrium configuration actually forms from generic 
initial data in BD theory. In this section, we answer this question by demonstrating 
the formation of boson stars in BD theory. The formation of boson stars in GR 
has been discussed by Seidel and Suen |^, for some set of initial data of non-zero 



measure. 

We start our evolution with initial data given by a Gaussian packet in the bosonic 
field $. This represents a local accumulation of matter field: 

$ = aexp(-6x^), (296) 

where a and b are free parameters. We set the BD field ip to be flat at the initial 
stage, so as to see if local inhomogeneity of the matter will form a boson star in BD 
theory. We reintegrate the lapse equation and Hamiltonian constraint equation to 
provide metric functons on the initial slice. We then observe its evolution. 

We find that, with particular parameters a and b, this system forms a stable 
equilibrium configuration, which can be recognized as the formation of a boson 
star. As a demonstration, we here show an evolution with parameters a = 0.1 and 
b = 0.025. The BD parameter uibd is taken to be 600. In Fig. 36, we show the 
BD field </) as a function of radial coordinate at various earlier times of evolution. 
We see that the BD field becomes negative quickly and begins oscillating around a 
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particular value. In the figure one can see that the BD field if between r = 40 and 
r = 50, and at around t = 50 is still not near its final configuration. The long time 
evolution of BD field is shown in Fig. 37, in which we show the BD scalar field ^bd 
at the center for early times, intermediate times, and late times. We can see the 
field settling down to a periodic oscillation in final phase, similar to the migration 
and transition cases in the previous section. The initial mass of this configuration 
in units of G^ is 0.39 and the final mass (at the end of our simulation) about 0.384. 
At this stage the magnitude of the central boson field is oscillating between 0.032 
and 0.048. The BD field oscillates between -0.00126 to -0.00166. This range of 
boson oscillations corresponds to masses between 0.342 and 0.410 respectively while 
the BD field oscillations give a mass between 0.355 and 0.405 respectively. Given 
that the mass at this stage is 0.384 (consistent with the above) we expect that the 
final mass will be between 0.355 and 0.384. 

We show the luminosity L{—dM/dt) versus time t curve in Fig. 38 during the 
same periods as in Fig. 37. In the early stage, we see one pulse is emitted from the 
system. This is related to the outgoing pulse from our initial boson field configu- 
ration. After this initial pulse, the system starts oscillating with a period tt (twice 
the frequency of the boson field oscillation) and the star begins forming. After that, 
we see the luminosity L begin damped oscillations as a function of t. (We cut out 
the initial large amplitude luminosity around t = 200). Even though the system's 
evolution is followed for a long time, the accuracy of the calculation is quite good 
with the Hamiltonian constraint satisfied to order 10~^ or better. 

We also note that certain parameters a and b (mentioned at the beginning of 
this section) will result in star formation but others will not. If we choose large 
amplitude a and small b, then the initial configuration has too large a mass and 
is not dispersive enough resulting in black hole formation during evolution. In the 
opposite limit if we have a very narrow localized wave packet it has a tendency 
to be dispersive (as per the wave equation). So if b is too large no boson star 
forms. Intermediate between these two are configurations that form stable stars. 
For example, if a = 0.1, then at 6 = 0.01 a black hole is formed, while at 6 = 0.025 
a boson star forms as we have shown, and at 6 = 0.035 the configuration disperses 
away to fiat space at the end of the evolution. On the other hand for b = 0.01 and 
a = 0.05 a boson star forms. This demonstrates that the boson star is a realizable 
object even in the BD theory, and opens windows to study them and other similar 
nontopological solitonic objects in astrophysical roles. 
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Quasinormal mode oscillation of a stable hoson sfar. The niaxiniuni value of the met- 
ric Qrr, and the central Brans-Dicke field ip{r = 0) are plotted as functions of time. Both 
of them take on the QNM frequency of the star. The oscillation is virtually undamped for 
a long period of time. 



Ill 




Figure 27 The magnification of the second figure of Fig. 26 is shown. In the diniensionless 
units we have adopted, the time interval between the ticks on the horizontal axis is it. We 
see the Brans-Dicke Held oscillates with period vr. Thus it has twice the frequency of the 
boson field, which oscillates with period 2 vr. 
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Figure 28 Finite perturbation of an S -branch boson star. Maximum metric grr is plotted. 
The metric is damped in time as the star settles to a new configuration. 
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Figure 29 The same model as Fig. 28. The central Brans-Dicke field (p{r = 0) is plotted. 
We see the oscillations damp out as the star settles down, indicating again transition to a 
new stable boson star configuration. 
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Figure 30 The same model as the previous figures. 'Luminosity' L = —dM/dt is plotted 
versus time. Clearly the radiation is decreasing in time. 
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Figure 31 Migration of an unstable boson star to a stable configuration; central Brans- 
Dicke field f{r = 0) and maximum of Qrr is plotted. There is a sharp initial drop in the 
radial metric as the star moves to the stable branch. The oscillations damp out in time 

as the star settles down. 
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Figure 32 Comparisons of total mass of the system M rescaled by its initial mass Minu 
during a migration process from U-branch star. Three lines are plotted. Although the 
mass-loss is similar in the three cases, the curve corresponding to the higher Brans-Dicke 
parameter is clearly closer to GR (as it should be). 
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Figure 33 Dynamical transition from U-hranch star to a black hole. The metric goo is 
plotted. The collapse of the lapse function is indicative of imminent black hole formation. 
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Figure 34 Dynamical transition from an excited state to a ground state boson star 
configuration. The metric Qrr is plotted. The initial four peaks indicative of a three-node 
star cascades to the ground-state branch after a long time evolution. 
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Figure 35 Dynamical transition from an excited state to a ground state boson star 
configuration. The metric grr is plotted at later times to show its oscillations after the 
star reaches a ground state. 
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Figure 36 An example of formation of boson star in Brans-Dicke theory. Snapshots of 
^BD{r) are plotted for initial stage of evolution. 
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Figure 37 An example of formation of boson star in Brans-Dicke theory. Dynamical 
behavior of the Brans-Dicke scalar field (Pbd{x = 0) is plotted for three evolution regions: 
early time, intermediate time, and late time. We can see the field settling down to an 
equilibrium configuration (periodic oscillation). 
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Figure 38 The emitted luminosity L = —dM/dt is plotted as a function of time. In the 
first stage of evolution the luminosity data takes on an underlying oscillation of period vr. 
The frequency of this oscillation is twice the frequency of the underlying boson field. This 
occurs after the emission of one Gaussian-shaped scalar pulse, related to the initial field 
configuration. The amount of mass loss decreases in time as the star forms and settles 
down. 
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3.5 Boson Halos 

In the last chapter we described the nature of the equihbrium configurations of Boson 
Halos, made of massless scalar particles. To study the formation of these objects 
we use the evolution equations {BD evolution equations with a = m = 0) and use 
the same code described in the first section of this chapter. Earlier studies, both 
numerical and analytical, have confirmed that self-gravitating objects with massless 



scalar fields cannot be compact ||6T|, ^. Thus taking a Gaussian type localized 
distribution of scalar matter and evolving the system using the evolution equations 
described above, results in dispersal of the scalar matter without forming any self- 
gravitating object. Even sinusoidal functions that were damped by exponential 
decays met with the same fate. These were functions like exp(— r) sin(r)/r and 
exp(— r) cos(r). In Fig. 39, a plot of the density versus the radius for different times 
is shown for a = 0.001 cos(r) exp(— r). The star dissipates very quickly. A Gaussian 
distribution represents a local accumulation of matter and is obviously the most 
likely happening. It seems though that only specialized distributions can give rise 
to these halos. 

The initial configurations for the scalar field that yield stable configurations were 
those characterized by a 1/r times a sinusoidal dependence at large r as well as an 
energy density p that had points of infiection. These were typically functions like 

^. (297) 

1 + r 

cos(r) 1 — exp ( — j , and 

sin(r) (l + -) log [l + ^— 
V r/ L 1 + r 

We show the density and radial metric evolutions for a field of the form 
0.003 cos(r)(l — exp(— 1/r)). This settles into a self-gravitating object after some 
time. Fig. 40a shows the density as a function of r. The central density p(0) in- 
creases during the evolution from its initial value at t = 0. Fig. 40b reveals the 
radial metric as it evolves in time as a function of radius. This too displays the 
settling to a stable configuration. In Fig. 40c we show the mass loss for the system 
as it finds itself a configuration. The mass at fixed radial values for these times is 
shown in Table 2. The amount of radiation for this system is relatively small and 
decreases in time. 

On the other hand, functions like l/(r + 1) , cos(r)/(l + r) and sin{r^'^^)/r^'^^) 
failed to settle to a bound state and just dispersed away. Scalar fields represented 
by these functions did not result in p having points of infiection. 
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3.5.1 Analytical Proof of Newtonian Stability 

To investigate the stability of Newtonian solutions under small perturbations using 
the perturbative method |^, ^, 64], we write the linear scalar field perturbation in 



the form 

$(r, t) := a(r)e-'"* + (5a(r)e'^"* , (298) 

where kn are the frequencies of the quasinormal modes. The perturbed scalar field 
equation 

6a" + 26a' /x + kl6a = , (299) 

together with the boundary conditions 

6a'{0) = 6a{R) = 0, (300) 

(for regularity of 6a at the origin and at the radius R of the boson halo) defines a 
Sturm-Liouville eigenvalue problem with real eigenvalues 

kl<kl< ... . (301) 

Eigenf unctions of this particular differential equation (299) are 

6a = ^^^^iM , (302) 

X 

where the eigenvalues are real, 

TT 

kn = n— ; (303) 

n is an integer and R the radius of the solution. Thus the eigenvalues k'^ are positive. 
This means that all modes are stable. 

3.5.2 Numerical Study of Stability 

We study the evolution of dilute configurations using the same numerical code we 
used for studying boson stars. We show a perturbed Newtonian configuration as 
it settles to a new Newtonian configuration. The perturbation that is used in this 
case mimics an annihilation of particles. A Gaussian bump of field is removed from 
a part of the star near the origin. Fig. 41a is a plot of the unperturbed versus the 
perturbed density at t = 0. The perturbed configuration is evolved and settles to a 
new configuration. The scalar radiation moves out as shown in Fig. 41b. In Fig. 41c 
the density profile is shown after the system settles down. The system is very clearly 
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in a new stable configuration. In Fig. 41d the mass is plotted as a function of radius 
for various times. Again one can see that the mass loss is decreasing by the end of 
the run showing that the system is settling down to a new configuration. The mass 
as a function of time is presented in Table 3 for different radii. 

We have so far only been successful in evolution studies of dilute configurations. 
The reason for this is that denser configurations need much better resolution in order 
for their evolution to be studied in our code. This is coupled with the difficulty that 
we still need the boundary to be very far away, so that the density has significantly 
fallen off, for an outgoing boundary condition to work. One possibility is to match 
our configurations with some other kind of matter at the boundary. For example, 
we could have a massless scalar field configuration surrounded by a massive field 
configuration like a boson star. Unfortunately, while we can find such configurations 
with smooth mappings of the metric and field density, the scalar field itself does not 
match smoothly for such configurations. So far, we are using either an outgoing wave 
condition or exact boundary conditions where the latter one simulates a vacuum 
energy. Further investigation is needed before we can decide whether these non- 
Newtonian configurations are inherently unstable. 

3.5.3 Radially Oscillating Solutions: Massive Field Case 

A boson star consists of scalar massive particles, hence in the simplest model one 
has a potential U = m^|<l>p. Exponentially decreasing solutions exist for special 
eigenvalues of the scalar field uj < m, so that the star has a finite mass. In the case 
oi uj > m, oscillating scalar field solutions can be found for all values of uj. The 
energy density reveals minima and maxima as opposed to the points of infiection seen 
for the massless case. Fig. 42a shows a comparison of the density profile for the two 
cases. In Fig. 42b we show the mass profile and the particle number as functions of 
central density. The binding energy is always positive making the system unstable. 
Numerically we find the system to completely disperse (as it should), making these 
parameters unfavourable as halo models. 
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Figure 39 The inability of niassless scalar Held configurations to form compact self- 
gravitating objects as discussed by Cliristodoulou and others ^ , [6^/ is verified numerically 
An initial field configuration of the form a = 0.001 cos(r) exp(— r) is seen dispersing in the 
plot. The dispersal takes place very quickly. 



Table 2 Mass values at different radii and different times for the evolution in Fig. 40. 
One recognizes that the change in the mass at fixed radius decreases in time, showing that 
the system is settling down to a stable conEguration. 
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Figure 40 {a) Top left: The formation of a self-gravitating massless scalar field object 
is shown here for an initial field conEguration of the form 0.003 cos(r)(l — exp(— 1/r)). 
The density in this case increases from its initial value as the conhguration settles down. 
(b) Top right: The radial metric for this configuration is shown evolving to a stable final 
conhguration. (c) Bottom: The mass is plotted as a function of radius at different times. 
As the configuration evolves it loses less and less mass as it settles down. 
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Figure 41 (a) Top left: An exact Newtonian configuration is perturbed and evolved using 
the GR code. The initial unperturbed and perturbed densities are shown. A Gaussian 
profile of scalar field is removed from the equilibrium configuration, and the constraint 
equations are reintegrated to give new metrics before beginning evolution, (b) Top right: 
The outgoing scalar radiation can be seen in a plot of the density at early times, (c) 
Bottom left: The density profile is shown after it starts to settle down. Between times of 
300 and 1500 very little has changed in the density profile, (d) Bottom right: The mass 
is plotted as a function of radius for differen^'^imes. In order to enhance the features only 
a small part of the radial region is shown. The amount of mass loss is clearly lessening in 
time. 
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Figure 42 (a) Left: The density profiles of massive scalar field configurations with to > m 
is compared to that of a massless held. The former has maxima and minima while the 
latter has points of inflection. This structure may be related to stability issues, (b) 
Right: The mass and particle number are plotted against the central density for boson 
star oscillators. We find in all cases that the mass is always greater than the particle 
number explaining their instability against a collective dispersion. 



Table 3 Mass values at different radii and different times for the evolution in Fig. 41. 
One recognizes that the change in the mass at fixed radius decreases, hence it proves the 
settling to a stable configuration. 
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Chapter 4 

Boson Stars in 3D Numerical Relativity 

We now describe our work on Boson stars in 3-D. One of the major drawbacks 
in going from a ID code to a 3D code was the paucity of grid points. Computer 
memory and time constraints hmited grid resolution. There was no simple outgoing 
wave condition (the boundary was no longer one point) and there was no reasonable 
way to implement a sponge. A sponge, while it prevents reflections, also uses up grid 
points. Besides, the length of the sponge had to be of the order of the wavelength 
of the outgoing radiation. Gravitational waves have longer wavelengths than scalar 
waves and this further inhibited the use of a sponge in a 3D code. 

The 3D code that had been developed to deal with vacuum spacetimes was 
extended to include matter terms. We used our model to test the code by comparing 
with established results in ID, as well as to results of 3-D perturbation studies. Once 
stability was achieved, we then used the code itself to get physical results of Boson 
stars under nonspherical perturbations. The main advantage with Boson stars as the 
matter source was that there were ordinarily no spacetime singularities as in black 
hole evolutions, as well as none of the surface problems that people dealing with 
neutron star configurations have to routinely face. The Boson field exponentially 
decays to zero at spatial oo. Nevertheless, the small amount of matter in the outer 
regions that prevents surface problems could, for the purposes of wave extraction, 
be perfectly well regarded as a vacuum. The code used a leapfrog scheme with the 
extrinsic curvature (related to the time derivative of the metric) and time derivatives 
of the Boson field on half integral time steps (n + 1/2) while metric fields and 
Boson fields were on the integral timesteps (n). To prevent grid point to grid point 
instabilities due to mesh drifting we routinely used a small amount of diffusion in 
the center of the star, as described in the last section of the introduction of this 
thesis. 

In an appendix we describe how we scaled the variables in our 3D code. This 
appendix shows how to relate a unit of computer time to physical time in terms of 
the mass of the configurations. 
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4.1 Stress Energy Tensor, Energy Density, Momentum 
Density, Sources 

The metric components in the ADM formahsm are given in (154) and (155). The 
lapse is A^ and the shift components are A^*. In the introduction to the thesis we 
wrote down the matter Lagrangian for a Boson star in scalar-tensor theory and 
showed how to get the stress-energy tensor from it. In GR we get the equations 
with the parameter a set to zero. The stress-energy tensor is then 



T, 



1 



/lU 



^1,M^1,!^- O^M^[^l,a^l'"+"^ ^1 ] + (^1 ^^2)- 



(304) 



where the subscripts 1 and 2 refer to the two real components of the complex Boson 
field. Since we rescale quantities to dimensionless ones no factor oim finally appears 
in the Einstein equations. (An appendix attached to this thesis explains the scaling 
in the 3D code). Thus we can write 



1 



(305) 



g'" (dt^PiY + ^1. (^^^ ^1. + g'^y ^1, + (7^^ ^1 
+ i^iy {g'' i^iy + g""' iJi. + g'^ ^u) 

+ ^1, {g^^ ^1, + g^y i,,y + g"^^ ^1 J + ^1 

Here the subscripts a;, y, and z refer to a partial derivatives with respect to those 
variables. The form of the density is given by n^n'^T^y where n° = 1/A^ and rt = 

-N'/N. Thus 



_ 1 T. N' N'N^ 



(306) 



The scalar field equation is the gravitational Klein-Gordon equation: 

= OA-V^ = ((7"^V^,/3);„-^ (307) 

1 



+ ^9"''9ua,pt^-^ 



1 



-9 



{V^),pt'-^ 



-9 



iv^tn,a-i^. 
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where the factor g represents the determinant of the 4 metric. Thus 

, ^ , U{N^-N^N,)^g-^i;^ -V^ = 0. (308) 

These second order (in time) equations {a = (3 = t) can be conveniently converted 



to first order by defining a new variable vr = J{N'^ — N^ Ni) (^/7 (?** ip) or, since we 
have two fields 



This gives the evolution of vr to be of the form 



''- = - Jp (Vl^a) a = 1,2. (309) 



TT 



^{N^-N'^N,)^^ (310) 



N' 



+ d,[N'7r)+dtU{N^-N^N,)^ — d,^ 



4.2 Testing the Scalar Evolver 

We computed the derivatives of the fluxes 7*-'\/(A^^ — N'^Ni)^ using a second- 
order scheme. We also adopted a second-order scheme for the first and second-order 
derivatives diipa and djipa- (The subscript a runs from 1 to 2 representing the two 
real components of the boson field). We used a staggered grid (the origin was not one 
of the points). In order to be better able to resolve the star we observed evolutions 
on an octant rather than the full grid so our perturbations also had to retain this 
octant symmetry. The first grid point was actually staggered half a spacing step 
away from the origin in one direction and the second grid point was half a step in 
the other direction. 

The first stage of the analysis was to get initial data. Spherically symmetric 
equilibrium data was obtained from the ID code, described in the last chapter. 
This data was then interpolated onto a 3-D grid. This was done by identifying the 
positions r{i,j, k) on the 3D grid, finding out which point it was closest to on the 
ID grid, and using the ID data before that point, at that point and the point after 
that, to get the value at r{i,j, k) with a three-point interpolator. This 3D data for 
the boson fields and for the metrics was then converted from r, 6, coordinates to 
X, y, z using 
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3ij ~ a- Q ■ 9i'j': 



x^ + y'^ + z^ = r'^, (311) 

X = r sin 6 cos </>, 
y = r sin 6' sin (p 

z = r cos9, and 
di[df_ 
di dj 

where the i and j coordinates refer to x, y, z and the prime coordinates to r, 6, 0. 

To test the scalar evolver we then froze the metric evolution and just allowed 
scalar field evolutions. (This is physically reasonable since an equilibrium boson 
star is characterized by a static metric.) When, after thousands of time steps, the 
boson field continued oscillating with a period of 27r (since the boson fields have 
time dependence of the form cos(tcode) and sin(tcode) where tcode = ^Et), we were 
convinced that our scalar evolver was doing well. 

The boundary conditions used were reflection symmetry at the origin and flat 
boundary conditions for the fields at the outer boundary. We present a plot of the 
scalar evolver in Fig. 43. 

When calculating {fip')' (where / is the flux and ip is the boson field) in the 
evolution of vr (equation (310)), it was numerically more advantageous to first write 

(/^')' = /< + />', (312) 

and then to calculate the two terms on the RHS separately. This avoids the numer- 
ical instabilities typically generated by finite differencing a finite difference. 

4.3 Equilibrium Boson Star 

The first step in testing the efficacy of the code was to take the known results of 
the ID code and reproduce them in the full GR code. Once the scalar evolver was 
found to be producing accurate results we wanted to test the code with full metric 
evolutions. Initially, we used equilibrium data obtained from our ID spherical code. 
The equations for full evolution have been developed in the introductory part of 
this thesis (|1.3.5|) . The four gauge degrees of freedom were used to set the three shift 



components to zero and then a maximal slicing condition was used for the lapse. 
4.3.1 Maximal Slicing 



The condition for maximal slicing pGl p5| is TrK = 0. Thus, on replacing for R 



from the Hamiltonian constraint equation (170) (with K = Q), into (174), we get 

dTr K 

— — = = -A^l^i, + N [KijK'^ + 4 7r(p + Tr 5)] . (313) 
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Here 5* = SI where Sij = Tij. The reason for writing it in this form was due to 
the presence of second derivatives in the Ricci tensor. The code had originally been 
developed for black hole spacetimes. Late time peak development in the metric 
functions rendered computation of second derivatives difficult. The above form was 



therefore numerically more advantageous |^ . This was important in the boson star 
collapse problem as well. 

One of the advantages of the maximal slicing condition was that it provided a 
linear equation for the lapse. The disadvantage was that an elliptic equation had to 
be solved on each time slice. Asymptotic flatness required the lapse to go to some 
constant in the limit r ^ oo. In our code the value of this constant was influenced 
by the underlying frequency of the scalar field and the lapse value at oo was the 
value of that frequency (see chapter 2). 

The main advantage of maximal slicing was its singularity avoiding ability. One 
of the advantages of a good coordinate system is avoidance of physical singularities 
by slowing down of evolution in the spatial region near such a singularity. The 
maximal lapse collapsed, vanishing exponentially, near the singularity and froze the 
evolution there (near r = 0, the physical singularity for a Schwarzschild blackhole 
or for a Kerr metric), allowing the maximal slices to span the rest of the future 
development of the initial data. This was important in the boson star collapse 
problem. Coordinate volumes usually shrink in the neighborhood of singularities. 
This is avoided by maximal slicing since the TrK = condition translates in the 
zero shift case to the volume ~ 7 being constant. Another advantage of the maximal 
slicing condition is that asymptotically it turns to a natural radiation gauge in which 
the gravitational radiation has a wavelike propagation. 

4.3.2 The Elliptic Solver 

In equation (313) there is a term of the form ^\Z\N and, as we did for the Klein- 
Gordon case (307), we write the maximal slicing lapse equation: 

d,{^Y'd,N) = N^ [K'^K^j+A-nip+S]) (314) 

The code used a conjugate matrix solver developed at the NCSA to which we 
added the matter terms. The basic technique was to first finite difference the above 
equation to second order. Because we were not deahng with flat space there were 
mixed derivative operators didj with i ^ j as well. Thus the second order finite 
differencing did not just connect points {i + 1, j, fc), (i, j, fc), {i — 1, j, /c), (i, j + 1, fc), 
(i,j — l,fc), ihjjk + 1), and {i,j,k — 1), a seven-point stencil. It also connected 
12 more points (i + 1, j ± 1, k), [i — 1, j ± 1, k), {i + 1, j, k ± 1), {i — 1, j, k ± 1), 
(i, j + 1, A; ± 1), and (i, j — 1, fc ± 1) where i j and k are the x, y and z coordinate 
positions respectively. Thus we were dealing with a 19-point stencil. For ease of 
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visualization we show the 2D case in Fig. 44 and the 3D case in Fig. 45. (In 2D 
we have a five-point stencil when there are no mixed derivatives and a nine-point 
stencil when mixed derivatives are present.) 
This gave an algebraic equation of the form 

p,g,r=+l 

/ . Q,i,fc Ni+p,j+q,k+r = (315) 



p,q,r=—i 

where the sum was over non-repeating terms and had 19 coefficients C. It had to 
be solved at each grid zone except boundaries, 

The conjugate gradient method involves solving the problem by solving an anal- 
ogous minimization problem. Consider a set of difference equations 

A/ = b (316) 

where / is the vector of unknowns (here our lapse value at each grid zone) and 
A is an A''^ x A''^ square matrix (A^^ = number of grid zones) containing the finite 
differenced coefficients at each grid zone. The right hand side is a vector of source 
terms. 

Define the quadratic function 

y(/) = i/tA/-6t/. (317) 

Therefore 

dV 

= Af-b = r. (318) 

of 

Thus solving (316) is the same problem as minimizing the potential V. Note that 
A is the matrix 

which is symmetric. So to apply this method the finite difference equations have to 
be written in symmetric form. 

The minimization of the potential function is carried out by generating a succes- 
sion of search directions s, with improved minimizers fi at each stage. In a maximum 
of Ng iterations a solution is found. 

However since we have a 19-point and not an iV^-point stencil, a lot of the matrix 
elements will be zero. (The matrix A is a sparse matrix.) Since the finite difference 
equations connect points immediately to either side of a a given point, as one steps 
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Figure 44 a) Top: The five-point stencil in 2D in Bat space is shown. The absence of 
mixed second-order derivatives in the x and y directions ensures that no points along 
the diagonal (except the center) are required, h) Bottom: In curved space the second- 
order derivatives include mixed derivatives (■^^)- We now have a nine-point stencil that 
includes the corner points. 
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Figure 45 a) Top: An extension of the previous figure to 3D is shown. The first figure 
shows a seven-point stencil, for the second-order derivatives in flat space, b) Bottom: In 

o2 fi'2 fi'2 

curved space the second-order derivatives include mixed derivatives ( q^q , Q^dz ' ^^^ dxdz ' 
and their permutations). We now have a 19-point stencil. 
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up the grid from point to point, there is a definite structure to system. The matrix 
has a banded diagonal structure with each diagonal band corresponding to one of the 
19 finite differenced coefficients. The coefficients are stored as 19 three-dimensional 
arrays, meaning 19A''„ not N"^ total elements. This reduces the number of steps 



needed to minimize the function. For details see 51, Bq 



4.3.3 Results of Maximal Slicing for Equilibrium Boson Stars 

With the above machinery in place we ran the code with maximal slicing. Our grid 
resolution was dx = dy = dz = 0.35 and the grid dimension was 96 x 96 x 96. The 
outer boundary condition on the boson fields was fiat (equal to the neighbor) and 
that on the metrics was static. In Fig. 46 we show a plot of the radial metric as a 
function of radius at various times. Clearly, this metric was not the static metric 
expected of an equilibrium star. 

The reason for this was a lack of gauge control. This is clear from the metric 
function ggg which should really be r^ in a perfectly spherically symmetric run. For 
ease of visualization we plot gee/f"^ as a function of r at various times and observe 
its drift from unity. This is shown in Fig. 47. 

Attempts to regain gauge control are described in the next two subsections. 

4.3.4 Role of the Shift Vector 

As we saw in the description of ADM in the introduction, we can write t^ as a sum 
of two terms, one along the normal to the hypersurface and the other the shift along 
the hypersurface. So far we have been using normal coordinates with shift vector 
equal to zero. Although this was the simplest choice of shift, with this choice, we 
have squandered our ability to implement gauge control by making specific choices 
of shift. From the plots of g^o it is clear that there is coordinate motion {dr/dt) 
due to numerical error. Maybe setting the shift terms to zero as we have done left 
too much residual spatial coordinate freedom? Might a shift term control this? The 
question now is how to choose a shift. To start with, the obvious choice seemed to 
be one that satisfied 

^ = 0. (320) 

dt 

For simplicity we used a spherically symmetric form, with 
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This yields, using the -^ = condition, 

2NKf,(, 
N'- = lilil^. (322) 

dPdee 

Typically this shift was phased in over a period of time. That is, a fractional 
amount of the calculated value of the shift was used in the code. The fractional 
amount was gradually increased from to 1. To circumvent dealing with boundary 
problems we isolated the points along the diagonal N^{i,i,i) into a ID array. We 
interpolated points to fill in values of r along this ID array that corresponded 
to the 3D grid spacing (grid spacing along diagonal = v3 grid spacing along x 
direction if the x, y and z directions have the same number of points). This was 
then extrapolated onto the 3D grid. This meant we needed only one boundary 
condition instead of worrying about all the boundary points on a 3D grid. Under 
our assumption of A^^ = A^*^ = then, 

iV^ = N'- -, Ny = N''- and A^^ = N'' -. (323) 

ly /y ry 

For the boundary condition we phased out the shift over the last 5 points 
{N''{nx-A) = 0.8N'inx-5), N'-{nx-3) = 0.6N''{nx-5) and soon till A^'^(nx) = 0). 
Although the shift rose to a peak near the center of the star, where gge had drifted 
the most before, and did have the intended effect of lowering the metric drift, it 
had a wavelike structure outside the star. The wave moved out and on reaching 
the boundary invalidated the code. A scheme of redefining the shift in the region 
where it first crossed the axis was tried. We used an exponentially damped attach- 
ment of the form A^*" = Aexp{—br), with A and b determined at each time step 
by matching the actual shifts at points close to where it crossed the axis with this 
expression. This and other variations on the functional form of the exponentially 
phased-out shift, lent a little additional stability to the metric but the metric even- 
tually developed a peak at the patching region of the shift and invalidated the code. 
Nevertheless, gee was fairly well locked in place and so it was at least certain that 
the drift had been due to coordinates drifting. 

A kicking scheme was also tried. The shift was fully phased in using the condition 
that best suited us at that point. A shift of the form ae~^^ , matched just before where 
the first shift peak crossed zero, was attached from the matching point onwards. 
After fully phasing in this shift, one discarded this method and tried to keep gee 
locked at its value at full phase in, by taking a snapshot of it at that time and then 
measuring the deviation from that value at each subsequent time step. The new 
shift at a given time at any point was either increased or decreased from its phased 
in value depending on whether gee decreased or increased from its phased in value 
at that point. Since gee ~ r"^ therefore ^.gee ~ 2r. Thus coordinate drift could be 
controlled by a shift of the form: 
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N^{t) = N^(t,) + c^^ (324) 

^^^ + 2r dt ' 
where tp was the phase in time and c could be positive or negative depending on 
which direction the drift was heading. In principle \c\ was unity, but in practice it 
tended to be a less and one had to tune this parameter. This kicking scheme had 
previously been used to lock the horizon in the black hole codes [|67 |. 



Although we found definite coordinate control with these methods the code al- 
ways became unstable at the outer boundary. We present comparisons of the radial 
metric function in simulations with maximal slicing with and without shift in Fig. 48. 
In Fig. 49 we present a comparison of gee with and without a shift phased in. Clearly 
the code was more stable with the shift phased in. Unfortunately the code ran longer 
without a shift. 

Even hybrid schemes were tried, with kicking done only in the outer region and 
the inner shift peak, that worked so well to lock gee-, retained. This led to grid point 
to grid point instabilities that diffusion could not cure. 

4.3.5 The K-Driver 

With the shift schemes providing some control, but not to the degree that was 
needed, it was timely that at this time the iT-driver was discovered. We were able 
to modify the i^-driver for the boson star problem p^ . 



The whole point was that we needed long term gauge control, and we had been 
using coordinate conditions that had been tried in numerical relativity for the study 
of black hole evolutions. The coordinate conditions themselves were becoming unsta- 
ble due to the accumulation of numerical errors during long evolutions. Shift vectors, 
which shifted coordinate points so they did not move normal to the constant time 
hypersurfaces, and time slicing conditions determining the normal direction of mo- 
tion of grid points, were gauge freedoms satisfying equations that were themselves 
subject to numerical inconsistencies. 

Consider the maximal slicing condition i^ = ^ = providing the elliptic equa- 
tion for the lapse (313). Any perturbation (numerical roundoff) in the lapse that 
made K nonzero, could not be forced by the slicing condition to get K back to zero. 
The lapse might evolve to its preferred value but errors made during its perturbation 
at some time in the past could still create errors in the slicing condition. The initial 
perturbation errors in the lapse might have been caused by numerical inaccuracies 
like errors due to finite differencing of nonlinear equations, with low resolution. The 
condition 

dK , , 
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written with K = could evolve to something where K ^ 0. Instead one considered 
a "K-driver" slicing condition 

OK 

-g^ + cK = 0, (326) 

which without K set to zero was 

cK = N^'ii - N [Kij K'^ - 4 TT (p + Tr 5)] - A^^ K^i. (327) 

While (313) did not actively enforce the K = condition the above slicing condition 
drove K to zero when it was perturbed away from it. The idea being that K ~ e"'^* 
exponentially went to zero, in the long run, for a positive c. One now could control 
the stability of the code by fiddling with the value of c that worked best. 

Finally, we actually used the K-driver for the boson star with shift equal to zero 
to get our best results. We needed a parameter c = 2 near the center (typically 
about a third of the grid) and then we let it smoothly damp away (^ce^s(^^^/^^ with 
R = r{nx,ny,nz) for all values of r > R/3). This damping factor was needed to 
prevent boundary problems from arising in the code. Since the drift was typically 
observed where the star was concentrated, this condition worked very well. We 
present a comparison of the zero shift maximal slicing versus the zero shift K-driver 
in Fig. 50. 

4.4 Spherically Symmetric Perturbations of a Boson Star. 

Once the 3D code had been stabilized to correctly display the true evolutions of 
an equilibrium boson star, we went to the next step of our code tests. This was to 
take a boson star and observe how the code displayed its evolutions under spherical 
perturbations. A well tested 1-D code provided us with initial data. The initial 
data was time symmetric and the momentum constraints were satisfied by ensuring 
that the perturbation kept ipi = ip2 and '02 = —''Pi = for the two boson fields. We 
show an oscillation of the radial metric in figures 5 1 and 52 (the former showing the 
downward part of the oscillation and the latter the upward part of it) of a boson 
star of original mass 0.5916 Mp^/m which has been reduced by about 10% after 
perturbation. 

A successful code test we performed was in calculating the mass loss of a boson 
star under spherical perturbations and comparing it to the results we had in ID. In 
the 1-D code we calculated the mass of the star using the ADM mass. Just before 
the sponge region (see chapter 2) we calculated the mass of the star from the radial 
metric 



M = - 1 . (328) 
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We could do this because in this region the field was small enough for the spacetime 
to be considered to be Schwarzschild. This was obviously not a good formula for the 
3D code with static outer boundary conditions. We needed a more global measure 
of the mass so that dominant contributions from the interior of the star would not 
be affected by outer boundary errors. In order to calculate the mass loss we used 
the fact that T*^ represents the energy flux {ergs /cm'^ sec) in the radial direction. 
Thus integrating T*'' over time and elemental area provided a measure of the mass 
loss at a given time due to scalar radiation. We calculated the mass loss in the ID 
code by using the ADM mass formula as well as integrated T*'" and compared the 
results to the 3D case. We ran the ID code with the same resolution as the 3D code 
and put detectors at the same position. We found very good agreement in all cases 
with different detectors and different types of perturbations. We show two plots, 
one with a localized spherical perturbation and the other with a global spherical 
perturbation, in Fig. 53 and Fig. 54 respectively. The latter plot also demonstrates 
the effect of resolution on the matching of the energy of scalar radiation. 

This was indeed the key code test that showed that the code was performing 
well given the computational resources available. 

4.5 Non-Spherical Perturbations of a Boson Star 

Finally we were ready to study the behavior of boson stars under nonspherical per- 
turbations. In the Introduction to this thesis we have already stated how Einstein's 
equations under nonspherical perturbations can give rise to a wave equation. These 
gravitational waves, along with the scalar radiation due to our scalar fleld, are the 
two big signatures of the boson star. We have seen in the Introduction, that the 
Zerilli functions, Newman-Penrose Spin Coefficient \E'4 and the Bel-Robinson vec- 
tor can provide information about the nature of this radiation and the gravitational 
waveform. One concern was of course whether one could extract the full waveform 
since the outer boundary effects could cause inaccuracies after a fairly short time as 
indicated by the results of the last subsection. Since the waveform extractions had 
to take place in the outer region of the star, which could effectively be regarded as 
a vacuum region, reflections from the boundary would cause problems earlier than 
if our detectors were placed well within the star. Here was another advantage to 
dealing with a boson star. As discussed in the introduction of this thesis (|1.3.3|) , 
in |^3| it was actually shown that the normal modes of boson stars damp out on a 
short time scale. The short time scale allowed us to measure all of the gravitational 
wave signal that contributes significantly to the gravitational energy output. 

We used a grid with one octant, making sure our perturbations were axisym- 
metric. We used an equilibrium configuration determined from the ID code and 
imposed on it such a perturbation. These were typically weak perturbations. Now 
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the constraint equations had to be solved on the initial time slice in order to provide 
consistent fields and metric components. This was done using an initial value solver 
that used a relaxation method to solve the elliptic equation. 

• Basics of the Initial Solver 

We used the initial solver to provide new metric components by solving the 
Hamiltonian constraint equation (170) on the initial time slice. York's procedure 



was followed |36]. We specified the initial scalar field configuration and an initial 
guess for the three-metric. The idea was to solve for a conformal factor cf) so as to 
get the real metric from the guess through the relation 

lij = <P%j (329) 

where the trial metric was 7jj. Assuming time-symmetric data, Kij = 0, the Hamil- 
tonian constraint equation takes the form 

R = 167rp. (330) 

With the updated metrics (in terms of the guess and the conformal transformation) 
we could calculate the connection coefficients in terms of those corresponding to the 
guess and also get the Ricci scalar in terms of the guess 

R = R(f)-^ - 8(j)-^did'(f). (331) 

The field configuration was retained as initially specified. This, along with •jij = 
(p'^^ij, was used to calculate p. Then the Hamiltonian constraint R = IGvrp was 
tested to some tolerance. So the problem was one of determining the conformal 
factor (j) satisfying 

Rcj)-^ - 8(j)-^did'(j) = 167rp(0, TT, (f)%j). (332) 

This could be solved using the CMSTAB (conjugate gradient method) by treating 
0^ X the right hand side (RHS) of the above equation as a constant for each iteration. 
The version we used in the code uses a Jacobi type relaxation technique. This was 
to write an equation for the conformal factor in the form 

dt(j) = e{LHS-RHS), (333) 

where LHS and RHS were the left and right hand sides of the constraint (conformal) 
equation (332). This gave the n + 1th iteration in terms of the nth one 

SO that one could relax to the solution and stop when LHS= RHS (which kept 
getting updated) was satisfied within some tolerance. Typically edt was chosen to 
be a fraction of the square of the grid spacing. 

Our perturbations were small enough that the initial guess for the metric was 
close enough to the actual value that this was not extremely time consuming. 
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4.6 Results of Non-Spherical Perturbations 

Far from the source, \l/4 represents an outgoing wave. It is thus normal to a two 
sphere of constant radius and the energy this wave carries over the whole sphere per 
second must be its integrated intensity. The Energy calculations from the Newman- 
Penrose scalar are then based on the formula 



E 



NP 



— I dt 



dt'^^A 



r'dVt. 



(335) 



We first estimated the degree of asymmetry in the system when subjected to 
spherical perturbations. These asymmetries could be due to numerical errors in- 
duced by, for instance, always finite differencing second derivatives in a particular 
order. For example taking the x derivative before the y derivative whenever one 



needed to find 



dxdy 



of any function. In Fig. 55 we show a plot of the comparison 



of gravitational energy measured at a detector for the case of a star that is unper- 
turbed, the same star under small and large spherical perturbations, and under a 
small nonspherical perturbation. 

The energy was also calculated by other means, namely the Bel-Robinson vector 



dE 
H 



1 

4:71 



I dt' (±^J\p^x + pyy + p^-2|/2 j 



r^dVt, 



as well as from the Zerilli functions \I'^: 

dE_J_ f d'^z 
'dt ~ 32^ \~dr 



(336) 



(337) 



We show a perturbed stable star with unperturbed central field density 0.1414 
that had been given a very tiny non-spherical perturbation so that we could obtain 
its quasinormal modes. Fig. 56 shows three plots of the gravitational energy output 
as detected by three detectors placed at various points on the grid. The first plot 
shows a Bel-Robinson calculation, the second a Newman-Penrose calculation and 
the third a Zerilli calculation. While the former two have consistent results the 
Zerilli plot shows much more sensitivity to detector positioning. That is, boundary 
effects as well as how close one is to being in a vacuum, seem to affect the Zerilli 
calculation more than the others. In the case of the Bel-Robinson and Newman- 
Penrose calculation the outermost detector is most affected by the boundary errors, 
while the other two give curves that flatten out. In Fig. 57 we show three plots, 
each now an energy comparison by the three methods, at a particular detector. The 
Bel-Robinson and Newman-Penrose calculations are most sensitive to how close 
one is to being in a vacuum while the Zerilli is most affected by boundary problems. 
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Taking the middle detector as the best compromise between boundary effects and 
degree of vacuum, we see that the Bel-Robinson and Newman-Penrose calculations 
match up very well for the most part until the point where they level off, with about 
12% differences at levelling off positions. 

The scalar radiation still dominates the system's energy output as seen in Fig. 58 
where we have plotted a tenth of the scalar radiation and the total gravitational 
energy radiation from Bel-Robinson calculations so as to be able to show them on the 
same plot. The matter radiation clearly moves out slower. The star's configuration 
after it gets rid of the gravitational radiation becomes spherically symmetric, and it 
continues its evolution through spherical oscillations. 

The nonradial quasinormal modes for the configuration with central field density 
0.1414 have been calculated in |]33| using perturbation theory. We have briefly 



described this in the introduction of this thesis. We have looked at our Newman- 
Penrose \l/4 function and performed a two-mode fit of their lowest two / = 2 modes. 
Our results are shown in fig. 59. Clearly we have a very close agreement to their 
results. After the first burst of energy, the star moves into the quasinormal mode 
and stays there till boundary effects become significant. 

We also compared two evolutions of a stable star, one with a small and the 
other with a large nonspherical perturbation. We observed that the energies of the 
gravitational waves and scalar radiation were greater in the case of the large pertur- 
bation. An infinitesimal perturbation caused the star to oscillate at the quasinormal 
mode of the original stable star while a substantial perturbation caused it to move 
to a new configuration, going through the modes of intermediate states. A Fourier 
transform for a star of central field density 0.15 with a small and large perturbation 
exhibits the differences in frequency in Fig. 60. 

Thus we have demonstrated that our code gives very consistent results and is 
capable of providing a detailed description of the gravitational signatures of a boson 
star. 

4.7 Colliding Two Boson Stars 

Another project, although still in its infancy, is that of the collision of two boson 
stars. The ID data of a single star was put on a 3D grid using reflection symmetry 
along an axis to simulate the presence of two identical boson stars. In order to center 
two stars on the Z axis separated by some number of grid points 2 nzs the boson 
fleld ip{i,j, nzs — k + 1) was added to ip{i,j, nzs + k) to get the new fleld ip{i,j, k) 
for all z, j, k on the grid. As an initial guess metric we just used the one star metric 
grrihjyk) = grrih j , nzs — k + 1) a.nd goo{i,j,k) = goo{i,j,nzs — k + l). We then used 
the initial solver to get proper initial data assuming the momentum constraints were 
satisfied (we had time-symmetric data and the n fields were defined appropriately). 
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The main result that we have right now, is that because the system produces much 
more gravitational radiation, the errors in the Zerilli functions due to boundary 
problems is much less significant. Typically, though, the first two detectors were in 
close agreement for the Zerilli function and the last two for the Newman-Penrose 
calculation. This shows the Zerilli function is still somewhat sensitive to outer 
boundary errors. The disagreement between the the two measures is markedly 
reduced compared to the earlier cases with smaller signals. So far we have only 
observed black hole formations from these collisions. These black hole formations 
were accompanied by negligible amounts of scalar radiation. Steep gradients in the 
metrics invalidated the code as they could not be resolved as the black hole formed. 
In Fig. 61 we show the gravitational energy radiated in the black hole formation 
during a two boson star collision. The central field density is 0.075 for the individual 
stars. One reason to study the two boson star collision is that it is a way to study the 
generalized two body problem. The amount of energy radiated may not sensitively 
depend on the component constituents of the collision. An order-of-magnitude cal- 
culation of the energy output to mass ratio (of the two boson star configuration) 
gives a value ~ 5 x 10~^ comparing well to the energy emitted (0.001) by two black 
hole collisions for black holes separated by large distances 

4.8 Boson Star: Black Hole Formation 

In the spherical collapse problem our ID code with the polar slicing condition caused 
the lapse to collapse and radial metric to grow uncontrollably at the approach of an 
apparent horizon. Our code was unable to handle these sharp gradients and crashed 
(uncontrollable numerical instabilities arose). By using a 3D code with maximal 
slicing, that was still singularity avoiding, but not as drastically as polar slicing, one 
could confirm the formation of the black hole. We took advantage of the ID code's 
ability to resolve the star, by monitoring the first part of its evolution in the ID 
code. Then after collapse was initiated, and before the gradients got too large, we 
collected data and placed it on the 3D grid. At this stage the star was smaller and 
the outer part of the grid could be conveniently thrown away allowing us decent 
resolution in 3D. 

We used "polar slicing" in the ID code, which did not set TrK = except on 
the initial time slice. As the system evolved in ID it evolved to TrK ^ 0. While 
maximal slicing [TrK = 0) avoided problems like grid stretching, by keeping the 
volume fixed in the absence of a shift, our initial data did not have TrK = 0. 
Although the K-driver transformed the slicing to i^ = after a short time the 
evolution exhibited severe grid stretching and we could evolve with zero shift only 
to a time of 22. 6M as opposed to t ~ 50M for vacuum Schwarzschild black holes. 
We present a plot of grr for the simulation of boson star collapse to a black hole in 
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Fig. 62. 

Singularity avoiding slicings offer a way to avoid evolution in regions with singu- 
larities while continuing the evolution in regions outside. However they only manage 
to delay the breakdown of a numerical code. Pathological properties of the slicing, 
like sharp gradients near the horizon, eventually invalidate the code. 

Since the region inside the black hole horizon cannot causally affect the region 
outside, one can in principle cut away the interior singular region and continue 
the evolution outside. An outside observer loses no information by cutting away a 
causally disconnected region. This requires an apparent horizon boundary condition 
(AHBC). 

For a collapsing star this can be implemented by evolving the star under suitable 
gauge conditions for a while and then introducing a shift vector after an apparent 
horizon is generated. This shift vector's function would be to maintain the apparent 
horizon at some fixed position ETf. This machinery had been developed for black 



holes [^ and was used in our code to lock the horizon in the spherical collapse 
problem. The simulation began without AHBC implemented and covered the r = 
region. After some evolution the apparent horizon formed and the grid was cut. A 
B-locking shift as previously described was applied to prevent coordinate drift. The 
lapse function was updated by a gamma-driver slicing. The details of these con- 
ditions and the causal differencing (using one-sided derivatives inside the apparent 
horizon) can be found in [^ and |^. With this machinery our code ran twice as 



long. We show the radial metric evolution in Fig. 
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Figure 46 Radial metric drift, with maximal slicing, for an equilibrium boson star of 
central field density 0.05. 
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Figure 47 The drift ofggg/r'^ from unity is shown, as an indicator of loss of gauge control, 
for a spherically symmetric equilibrium hoson star configuration. 
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Figure 48 The radial metric function, showing markedly more drift at the same time 
under conditions of no shift, compared to a gee freezing shift with a kicking scheme. 
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Figure 49 A plot of ggg shows markedly more drift at the same time under conditions of 
no shift compared to a shift with a kicking scheme. Both are with maximal slicing. 
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Figure 50 A comparison of the metric function drift in the time evolution of an equi- 
Ubriuni boson star with maximal slicing and K-driving condition is shown. The latter is 
clearly capable of controlling metric drift. 
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Spherically perturbed Boson Star 

Downward Part of Radial Metric Oscillation 
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Figure 51 We show the downward part of the radial metric oscillation of a boson star 
that has been perturbed by a spherically symmetric perturbation. As the star expands 
the radial metric descends. 
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Figure 52 As the star starts to contract after its expansion in reaction to the perturbation, 
its radial metric starts to rise again. This part of its motion is shown in the figure. At the 
last step we start to see it reverse its motion again. The central field density of this star is 
■01(0) = 0.15 in the dimensionless units discussed in chapter 2. Its mass is 0.5916 Mp^/m. 
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Figure 53 The mass loss in ID computed using the ADM mass formula, and integrated 
energy Hux T**", is plotted along with the integrated energy flux T**" using the 3D code 
for a local spherical perturbation. The closeness of values is remarkable until the effects 
of boundary reflections starts affecting the accuracy of the 3D code. 
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Figure 54 The mass loss in ID computed using the ADM mass formula, and integrated 
energy Eux T**", is plotted along with the integrated energy Eux T**" using the 3D code for 
a global spherical perturbation. The closeness of values is remarkable until the effects of 
boundary reflections starts affecting the accuracy of the 3D code. In addition the effect 
of grid resolution is shown by also plotting a lower resolution 3D calculation which shows 
greater discrepancy from the ID result but which is still quite good. 
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Figure 55 A comparison of gravitational energy output for a star of unperturbed central 
field density 0.15 for spherical and nonspherical perturbations. While the level of this 
'background' is not insignificant it is sufficiently small that the results for the nonspherical 
perturbations may be trusted. 
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Figure 56 The gravitational wave energy output of a non-spherically perturbed stable 
boson star of central field density 0.1414 is shown as measured at three detectors at the 
outer region of the grid. Three plots are shown, the first showing measurements using the 
Bell-Robinson vector, the second using the Newman-Penrose ^4 parameter and the third 
using the ZeriUi functions. The most erratic are the ZeriUi measurements. 
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Figure 57 A comparison of the gravitational energy output using the three measures 
defined in the previous figure is made at three detectors for the non-spherically perturbed 
stable star described in the previous figure. While the Bel-Robinson and Newman-Penrose 
^4 calculations are quite close the Zerilli data is more erratic. 
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Figure 58 The stable star previously described, when subjected to a nonspherical pertur- 
bation, emits gravitational and scalar radiation. The gravitational radiation is far less in 
amount compared to the dominant scalar energy radiated. The scalar radiation, due to the 
bosonic matter, travels slower than the gravitational wave due to curvature perturbations. 
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Neuman-Penrose Scalar ^^ 

Nonspherical Perturbation of Stable Boson Star 




Figure 59 The gravitational waveform ^4 is shown, as is a quasinormal mode fit for 
the stable hoson star previously described. The star gradually starts oscillating in its 
quasinormal mode and then at late times boundary effects start affecting the system. The 
plot has time in units of the underlying oscillation of the scalar field (time period of2ir). 
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Figure 60 The nonradial modes of a strongly perturbed star are shown to be very different 
from the nonradial modes of slightly perturbed one. The former star will lose far more 
mass through scalar and gravitational radiation and move to a new configuration. The 
latter oscillates in the quasinormal mode of the original configuration. 
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Figure 61 The energy carried by a gravitational wave as measured by three detectors, is 
shown for a Zerilh calculation and a Newman-Penrose calculation. The former shows very 
close agreement for the Erst two detectors but outer boundary effects result in inaccuracies 
in the third detector's signal. The latter on the other hand shows a close agreement in 
the last two detectors' signals. 
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Figure 62 Tin's piot siiows tiie evolution of the radial metric function §„ for the case of 
the spherical collapse of a boson star to a black hole. In this simulation the shift is zero 
and the lapse function is obtained by solving the maximal slicing equation with a K-driver 
term ~ cK. 
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Figure 1: The evolution of the radial metric function grr is shown for the case 
of a perturbed boson star collapsing to form a black hole. The formation of 
the apparent horizon occurs around t = 5.6M. The evolution of the resulting 
black-hole spacetime is continued up to t = 45. 3M through the implementation 
of an apparent horizon boundary condition. The evolution in the inner region 
is quite stable with the observed drifting attributable to the error in solving the 
elliptic equation for the lapse function and errors from the outer boundary. 
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Chapter 5 
Conclusions 

Even the skeptic, who refuses to beheve in the possibihty of the existence of self- 
gravitating bosons, would be forced to admit that they are indeed very useful objects: 
That is to say, they are rich in physical content, providing insights into the nature 
of several physical systems. 

For instance the dark matter problem is real and there must be something to 
explain it. We know that normal baryon densities are not enough to close the 
universe, and we know that cold and hot dark matter models fail to explain large 
and small scale structure of the universe, respectively. Relativistic Bosons seem to 
offer a partial solution to these problems. Although they have particles with high 
velocities to explain the large scale power spectrum, they have more low velocity 
particles than fermionic hot dark matter (because of the Bose-Einstein distribution 
as opposed to the Fermi-Dirac distribution). Thus while small scale power is erased 
for hot dark matter due to free streaming, enough survives on the smallest scales 



with this kind of bosonic hot dark matter [Q. 

Also, many particle physics and cosmological models rely on some form of scalar 
particles. These unseen scalars could very well be dark matter constituents. From 
the primordial soup of scalar particles self-gravitating objects can form by the Jeans 
instability mechanism. Perturbations larger than the Jeans' length on a homoge- 



neous background could cause the compactification of matter [|T6 . 

Self-gravitating objects made from complex scalar fields are called Boson stars. 
These solutions to the Klein-Gordon-Einstein equations are held together by the 
balance between attractive forces of gravity and the dispersive forces of the wave 
equation. In the literature they are referred to as mini-soliton stars because of 
their small sizes. The addition of a small coupling parameter can however increase 
their size to neutron star dimensions, making them astrophysically interesting. 

The masses of these stars increases with central density up to a peak and then 
decreases. We have investigated the stability of ground and excited state boson 



stars with self-coupling [Q. The branch to the left of the peak is called an S or 
stable branch while the branch to the right is called the ?7-branch. Configurations 
on the S branch have very specific quasinormal modes of oscillation under small 
perturbations. Under finite perturbations they settle down to new configurations 
on the same branch. When perturbed, configurations that lie on the unstable or 
[/-branch either migrate to new configurations on the ^-branch or collapse to black 
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holes. These characteristics are shared by boson stars with or without self-couphng. 

Excited states are configurations with nodes. The field of an n^^ excited state 
star has n nodes and its radial metric has n+ 1 peaks. Their mass profiles are similar 
to the profiles of boson stars in the ground state, which makes it appear as if they 
have a stable and an unstable branch of configurations. However, irrespective of 
which branch they lie on, excited boson stars are unstable with different instability 
timescales. Low density excited stars having masses close to ground state config- 
urations will form ground state boson stars after evolution. Denser configurations 
form black holes with the decay time decreasing with increasing central density till 
one approaches the density corresponding to zero binding energy. As the central 
density approaches this critical central density the kinetic energy of the star starts 
to increase as it becomes more dispersive. It still collapses to a black hole but on a 
larger time scale. Beyond this point, for densities corresponding to positive binding 
energy, the stars disperse to infinity. 

An interesting feature in the collapse of excited state boson stars is that, during 
this process, they cascade through intermediate states, rather like atoms transiting 
from excited states to the ground state, suggesting that boson stars behave in some 
ways like gravitational atoms [Q . However, for boson stars, an investigation of the 
possible decay channels (selection rules) seems much more difficult (if at all possible 
or meaningful) due to the intrinsic nonlinearity of the theory. 

We have also studied the stability and formation of boson stars in Brans-Dicke 
theory |5^. The arbitrariness of units of length, time, and masses, and the observed 



coincidences in the values of dimensionless constants, resulted in alternative theories 
of gravity, like BD theory and scalar-tensor theories, being propounded. While it has 
been established that the weak-field differences between these theories and gravity 
is insignificant at the present times, the question is whether a boson star model can 
be a source of strong-field differences. We have worked on a BD model, with the 
view of carrying on the calculation to more general ST theories in the future. 

We find the basic features of [/-branch and ^-branch stability and excited state 
transitions to be similar to GR. We do see the BD field, which is static for an 
equilibrium configuration, begin to oscillate when perturbed, as do the metric com- 
ponents. However these oscillations of the gravitational field are not accompanied 
by any significant gravitational radiation. 

Using a Gaussian initial boson scalar field configuration, to represent the local 
accumulation of matter, we were able to demonstrate the formation of boson stars 
in BD theory. This demonstrates that a boson star is a physically realizable object 
in BD theory. 

The fiatness of galactic rotation curves indicates the presence of dark halos ex- 
tending beyond the last visible point, with masses proportional to distance and den- 
sities inversely proportional to distance squared. While boson stars in the ground 
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state have too fast a mass fall off to fit galactic rotation curves, third and fourth 
excited state stars with and without self-coupling have been used to fit them quite 



well |^8[, 1^. These states are unstable but their energy differences from lower 
excited states are large enough to warrant large instability time scales. As we have 
seen in our own simulations, a repulsive self interaction term further increases this 
instability time scale. 

Another model that can explain galactic rotation curves involves massless scalar 



fields [|^, [^. The dispersive effects of the scalar field prevents collapse of these 
halos to black holes. These configurations exist for different values of the scalar 
field oscillation frequency (0 ~ e*'^*) and for each value of central density. While 
the value of the central density is related to the orbital velocity, the scalar field 
oscillation frequency determines the mass and energy density. The Newtonian con- 
figurations had already been described and we studied the stability and formation 
of these objects in a relativistic codeHSI. The mass and energy densities still had 



the desirable properties for flattening out rotation curves. 

We have also shown that massless complex scalar fields are able to settle down 
to a stable configuration. It seems that the formation process needs a special profile 
for the energy density. The appearance of points of inflexion within the decreasing 
density supports the 'birth' of the boson halos. In contrast, extremal values for the 
density cannot be compensated and lead to the destruction of the initial configu- 
ration. Stable configurations cannot be formed from initial Gaussian distributions. 
This could mean that the formation of boson halos could be difficult. However under 
perturbations a known stable solution settles to a new one. 

The thing that makes these models undesirable though is that they are not 
asymptotically fiat. The cut-off radius of these configurations is quite arbitrary. It 
could be that the size of the halo depends on the value of the cosmological constant. 
We are working on boundary conditions to make this model more viable. In the 
process of studying these solutions we have found similar solutions for the massive 
scalar field. These solutions have energy density characterized by a series of maxima 
and minima rather than saddle points. They have positive binding energy and as 
expected disperse in our numerical simulations. 

We have furthermore investigated boson stars in a general three-dimensional 
code demonstrating the key role that boson stars play in the field of numerical 
relativity. By exploiting the difference between boson stars and neutron stars, we 
have managed to deal with them numerically. There are no surface problems in a 
boson star system. Nevertheless the similarities in their structures has allowed us 
to study stellar dynamics using them as a model. 

On the numerical front we have developed a fully relativistic code with matter 
sources. By comparisons to known results, we have been convinced that the code 
is indeed stable within the computational resources available to us. We have tested 
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our scalar evolver by itself and seen that it is stable for thousands of timesteps. We 
have seen metric components remain static for equilibrium boson stars, for several 
time steps before outer boundary problems invalidate the code. The amount of 
scalar radiation under spherical perturbations matches very well with our ID code 
until boundary reflections become significant. We have observed one full metric 
oscillation of the quasinormal mode under spherical perturbations. 

We have then moved a step further to study real physics with this code by 
studying the behavior of boson stars under non-spherical perturbations. We have 
measured the gravitational wave output using different techniques with comparable 
results. We have calculated the quasinormal modes of the star under infinitesimal 
perturbations. Our results are consistent with those of ^|. The nonradial modes 



for boson star configurations are strongly damped. This allows us to get the full 
gravitational wave output for these stars in a relatively short evolution timescale. 
After the emission of gravitational waves the star becomes spherical and then os- 
cillates in its spherical quasinormal mode. The amount of scalar radiation is much 
more than the gravitational radiation for these configurations. The scalar radi- 
ation decouples from the gravitational radiation. While the latter propagates at 
the speed of light the scalar radiation is slower and arrives later at our detectors. 
Under larger non-spherical perturbations the star emits more radiation and does 
not oscillate in the quasinormal mode of the initial configurations. By comparing 
the amount of gravitational radiation for unperturbed, spherically perturbed, and 
non-radially perturbed configurations, we have convinced ourselves that numerically 
induced asymmetries are negligible. 

We have collided two boson stars and formed black holes. This is with the 
view of studying the generalized two-body problem in GR. The amount of scalar 
radiation in these black hole collapses is negligible and the gravitational radiation 
dominates the system. The ratio of energy radiated to initial mass of the system 
is of order 0.0005 as opposed to 0.001 for two black hole collisions initially largely 
seperated. These numbers are comparable and might suggest that the two body 
problem may be independent of the nature of the constituent objects. Thus a boson 
system might, with its absence of singularities and surface problems, be ideal to 
study such systems. 

We have also seen the collapse to black holes of unstable boson stars. This is of 
importance on two accounts. When we performed these simulations in ID we used a 
polar-slicing condition which caused the lapse to collapse and radial metric to blow 
up at the approach of an apparent horizon. Thus we could not carry the simulation 
through to actual black hole formation . By using a maximally-sliced 3D code we 
can carry our simulation further and confirm black hole formation. In addition by 
running the ID code with very fine resolution until gradients start getting steep 
and collapse has been initiated, we could throw away the outer part of the star (its 
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radius has been reduced due to collapse). It could then be placed on a 3D grid with 
finer resolution than if we started the simulation from the beginning in 3D. After 
this we could turn on the AHBC and evolve the resulting black hole with matter 
present in the space-time. This added a new dimension to the black hole problem. 

Thus with a simple system of bosons, we have touched upon numerous physical 
and numerical issues. We now have stable ID and 3D codes that are capable of han- 
dling different system configurations. There are several models whose equilibrium 
configurations have been described but whose evolutions have not been studied. 

For example we have considered boson stars which are minimally coupled to 
gravity. If one considers the relevance of boson stars in the early universe, it would 
be worth investigating scalar fields non-minimally coupled to gravity. Since the 
Pecci-Quinn symmetry |]70| is not exact, QCD effects place the coupling strength 
for axions at ^ ~ 10~^° (not 0). This concept might extend to complex scalar fields. 
In |7T| a nonminimal coupling of the form ^ |0pi? is chosen and the equilibrium 



configurations are described. For a coupling strength ^ > 4.0 the mass is always less 
than the particle number. One expects the branch to the right of the maximum mass 
to be unstable but if the particle number is greater than the mass this instability 
must always manifest itself in gravitational collapse and not dispersion. These things 
can only be confirmed in an evolution code. 

There are also descriptions of the equilibrium sequences of mixed fermion-boson 
stars W% 173]. Pure boson stars are qualitatively quite similar to pure fermion stars. 



A mixture of bosons and fermions interacting only by gravitational forces is of 
interest in field theory as a study of the self- confinement of two quantum fields [ [T^ . 
Also bosonic axion fields could be captured by neutron stars and infiuence their 
stability. Similarly nuclear matter in high density neutron cores could be subject 
to phase transitions into the pion condensate. Again the stabilities and behavior of 
the system can be studied only in an evolution code. 

One interesting project would be to confirm whether real scalar field configu- 
rations could exhibit long term stability under non-spherical perturbations. While 
these configurations had no static equilibrium solutions |]73|, |7^, 0, periodic solu- 



tions have found ^Jj- Evolution studies in a ID code showed that there were oscillat- 
ing configurations that oscillated about some configuration for a long time, without 
observable diminishing mass. This study could be taken to the next step by studying 
the behavior of the system under more realistic non-spherical perturbations. Since 
axions are described by real fields this would be of particular importance. 

The list is endless and one hopes that the tools described and developed in this 
thesis will make their contribution in furthering knowledge in physics and numerical 
relativity. 
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Appendix: Scaling 

In transferring initial ID data from the ID boson star code, to the 3-D "G" code, 
where it was interpolated onto a 3D grid and evolved, one had to make sure all 
coordinates were rescaled appropriately. Essentially this meant going from a code 
where h = 1 and c = 1 (the ID code) to a code where G = 1 and c = 1. In terms 
of coordinates this translates to the spatial coordinates going from units of inverse 
mass to units of mass. Consider for example the ID Hamiltonian constraint. 
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In dimensionless coordinates (r, t, a, and A^) 
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where m refers to the mass of a boson and uiq the oscillation frequency of the boson 
field, this becomes 
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Since the field is very small at the edge of the grid the metric components can be 
taken to be Schwarzschild there. So g^{oo) ~ 1/(1 — 2MG/r). This means the mass 
of the star is given by 



M = JL 1-1 



2G 
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(341) 



calculated at the edge of the grid. 

Using ^ = Mpi (for the ID code therefore G = 1/Mpi) and r = r/m, where r 
is the dimensionless radius read off from the output of the ID code, we get 
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Going back to equation (339), one sees that when the code runs for a certain di- 



N 



tcodc it corresponds to physical time t = l/c^o^codc = ^^i^code, 



tr 



mensionless time tc 

where uo has been written in terms of the lapse. For a boson star of mass 
0.633 Mpi/m (maximum mass of a ground state star) we can compute the phys- 
ical time 
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where we have written the mass of the boson in terms of the mass of the star 
(Schwarzschild radius of the star ~ 1/m). Computing ^ at the edge of the grid, 



where N ~ J\ — 2M/r, we find this fraction to be of order 1. Thus 1 unit of code 
time in terms of the mass of the star is about 1.6 GMje'. 

In the ?)D code however G=l but h, is not. So consider the Hamiltonian constraint 
equation again. There is no special reason for this but it suffices to use it to figure 
out the scahng parameters. Every derivative must be multiphed by /i, which had 
been set equal to 1 in the \D case. Now G is set equal to 1. So one has for the 
Hamiltonian constraint equation 
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substituting dimensionless coordinates 



mr/h, t = ujot, a = yATchcj), N = 'N 
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we get back equation (340). Now to interpret our results, which have thus far been 
stated as raw dimensionless numbers: Without metric evolution the code runs to a 
time of 2000. The more complete problem must include metric evolutions. For this 
case, when we have in our reports been stating that the code runs to a time of 140, 
we actually mean that the code runs to a time of 140/a;o with 1/ujq = hlnirj^. In 
G = 1 units we have 
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calculated at the edge of the grid. This gives 
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Using the fact that G 
before. That is 
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1 and c = 1 implies % = Mp^ one gets the same result as 
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which for the problem in consideration is M = 0.05 Mpi/m. Replacing for 1/m in 
the equation for I/uq gives I/cjq = n(oo) cT^ where the h cancels the factor Mpi. N 
is > 1 and N < 1 with the ratio slightly greater than 1. So an order of magnitude 
for 1 time step in terms of mass for the maximum mass boson star is 1/0.633 ~ 1.6 
times the mass of the star. 

For our equilibrium configurations with full evolutions we use star masses of 
order 0.5 Mpjirn so one time step is slightly more than two times the mass of a star. 
The code runs to a dimensionless time of over 500 or over thousand times the mass 
of the star (tphysicai > lOOOGM/c^). 

For non-spherical perturbations the mass was of order 0.584 Mpi/m. One time 
step was roughly 1.87 times the mass of the star. 
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